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Abstract. We study hydrodynamic instabilities during the first seconds of core collapse supernovae by means of two- 
dimensional (2D) simulations with approximative neutrino transport and boundary conditions that parameterize the effects 
of the contracting neutron star and allow us to obtain sufficiently strong neutrino heating and, hence, neutrino-driven explo- 
sions. Confirming more idealised studies as well as supernova simulations with spectral transport, we find that random seed 
perturbations can grow by hydrodynamic instabilities to a globally asymmetric mass distribution in the region between the 
nascent neutron star and the accretion shock, leading to a dominance of dipole (/ = 1) and quadrupole (/ = 2) modes in the ex- 
plosion ejecta, provided the onset of the supernova explosion is sufficiently slower than the growth time scale of the low-mode 
instability. By gravitational and hydrodynamic forces, the anisotropic mass distribution causes an acceleration of the nascent 
neutron star, which lasts for several seconds and can propel the neutron star to velocities of more than 1000 km s _I . Because the 
explosion anisotropies develop chaotically and change by small differences in the fluid flow, the magnitude of the kick varies 
stochastically. No systematic dependence of the average neutron star velocity on the explosion energy or the properties of the 
considered progenitors is found. Instead, the anisotropy of the mass ejection and thus the kick seems to increase when the 
nascent neutron star contracts more quickly, and thus low-mode instabilities can grow more rapidly. Our more than 70 models 
separate into two groups, one with high and the other with low neutron star velocities and accelerations after one second of 
post-bounce evolution, depending on whether the / = 1 mode is dominant in the ejecta or not. This leads to a bimodality of 
the distribution when the neutron star velocities are extrapolated to their terminal values. Establishing a link to the measured 
distribution of pulsar velocities, however, requires a much larger set of calculations and ultimately 3D modelling. 
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1. Introduction 

Spect ropolarimetry (Leonard et all I20061 IWang et all I2003I 
l200ll and references therein) indicates that global anisotropies 
are a common feature of many core-collapse supernovae (SNe). 
Since the asymmetry inferred from the observed polarisation 
grows with the depth one can look into the expanding su- 
pernova ejecta, the origin of the anisotropy seems to be in- 
trinsically linked to the mechanism of the explosion. Recent 
high-resolu tion imaging of SN 1987 A with the Hubble Space 
Telescope llWang et all 12002) as well as large values for the 
measured space velocities of young galactic pulsars may be in- 
terpreted as a support of such a link. The average pulsar veloci- 
ties are as high as 200-500 km/s, and some neutron stars (NSs) 
move through interstellar space with more than 1000 km/s (e.g. 
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Cordes et alJI 19931 Lvne & Lorimerlfl99l lHansen & Phinnevl 
1997t IZou et al.1 12005b IChatteriee et al.1 l2005t) . Claims of a 
bimodality of the pulsar velocity distribution are still con- 
troversial. While some authors have obtained evidence for 
such a bimodality JCordes & Chernoffll998UFrver et all 19981: 
lArzoumanian et alJ 120021 iBrisken et alJ l2QQ3|). others found 
that a simple Maxwelli an fit works best (lHansen & Phinnevl 
ll997HHobbs et all2005h . 

Binary disruption, e.g. as a consequence of the SN explo- 
sions which give birth to the NSs, does not lead to sufficiently 
high velocities. Furthermore, the orbital parameters of many 
binary systems imply an intrinsic acceleration mechanism of 
the pulsars, pr obably linked to their creation (see lLaill20oil 
lLai et alJl200ll for reviews). Quite a number of explanations 
have been suggested, mostly involving anisotropic mass ejec- 
tion in the SN explosion or anisotropic neutrino emission of 
the cooling, nascent NS. The former suggestion might be sup- 
ported by the fact that some pulsars seem to propagate in a 
direction opposite to mass distribution asymmetries of their as- 
sociated SN remnants. But a clear observational evidence is 
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missing, and hydrodynamic simulation s have previously eithe r 
produced rather small recoil velocities (Janka & Miiller 1994), 
or started from the assumption that a dipolar asymmetry was 
already present in the pre-collapse iron core and thus gave rise 
to a l arge anisotropy of the SN explosion ( burrows & Havesl 
1996j). The origi n of such big pre-col lapse perturbations, how- 
ever, is not clear (Murphy et al 12004 . 

Suggestions that a "neutrino rocket engine" boosts NSs 
to large velocities (e.g. IChugailll984 iDorofeev et aljfl985t 

pplT ail 



Burr ows «fe_J^^c^slgv||l_g^6t IWboslevll 1 987t for reviews see lLail 
20011 and Lai et alJkOOll) make use of the fact that the huge 



reservoir of gravitational binding energy released during the 
collapse of the stellar core is mostly carried away by the neutri- 
nos. Creating a global emission anisotropy of these neutrinos of 
even only 1 % - which is sufficient to obtain a NS recoil of about 
300 km s -1 -, however, turned out to be very difficult. Most 

ideas refer to unknown neutrino properties (e.g., Fuller et alJ 

— ^ I II — " — | 

2003; Frver & Kusenko 2006 and refs. therein) and/or require 

the presence of a very strong magnetic field with a large dipole 
component (instead of being randomly s tructured and vari able 
with time) in the ne wly formed NS (e.g. lArras & La i 1999a b; 
ISocrates et al.ll2005l) . Such assumptions are not generally ac- 
cepted and are not the result of self-consistent calculations but 
put into the models "by hand". 

If the observed high pulsar velocities indeed go back to 
the early moments of the SN explosion, the simplest explana- 
tion would certainly be a common origin of explosion asym- 
metries and pulsar acceleration. In this case anisotropic ejec- 
tion of mass would lead to a recoil (or "kick") of the NS due 
to (linear) momentum conservation. Various kinds of hydrody- 
namic instabilities might in fact be responsible for a large-scale 
deformation of the ejecta and globally aspherical explosions. 
Perturbation analy sis of volume- filling thermal convection in 
a fluid sphere by Chandrasekhar (1961), found largest growth 
rates for the I = 1, m = mode (in terms of an expansion in 
spherical harmonics Yf of order /, m). This is supported by 
(full 4jt) three-dimensional simulations of c onvection in red 
giant and non-rotating m ain sequence stars dWoodward et alJ 
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2003|lKuhlen et alJ2003T) . Motivated by Chandrasekhajs anal- 
ysis, Herant ( 1995) speculated about the formation of a stable 
/ = 1, m = convective mode in the neutrino-heated layers 
between the gain radius and the supernova shock. In this con- 
figuration there exists only a single buoyant bubble (outflow), 
and a single accretion funnel (inflow), which reaches from the 
postshock region down to the NS. Herant ( 1995) suggested the 
potential importance of such a convective pattern for NS kicks 
up to nearly lOOOkm/s. Instability of the accretion shock to a 
global Rayleigh-Taylor mode which could lead to asymmetric 
shock expansion and a net recoil of the N S of several lOOkm/s 
was also predicted bv lThompsorJfc oOO). However, according 
to the linear analysis by Foglizzo et al. (2006), advection tends 
to stabilise the growth of long-wavelength perturbations in the 
neutrino-heated accretion flow behind the standing shock. A 
convective trigger of such instabilities therefore requires the lo- 
cal growth rate to exceed a critical threshold value. 

Non-radial, low-mode instability of shocked accretion 
flows can also be caused by the "advective-acoustic cycle". 
In the astrophysical context, this instability was first discussed 



by iFoglizzo & Tagged feOOOh and iFoglizzol i l200ll l2002h in 
application to Bondi-Hoyle accretion of black holes. It was 
more recently considered f or supernova-core like conditions by 
Idalletti & Foglizzol d2005h . This instability relies on the fact 
that the infall of entropy and vorticity perturbations produces 
acoustic waves that propagate outward and create new entropy 
and vorticity perturbations when reaching the shock, thus clos- 
ing an amplifying feedback cycle which eventually results in a 
dominant / = 1 or / = 2 mode. The advective-acoustic cycle 
can even operate under conditions, in which convective insta- 
bilities are hampered, e.g. if the advection of matter out of the 
convectively unstable region is too fast to allow for a significant 
co nvective growth of sm all perturbations. 

iBlondin et alJ J2003I) investigated numerically an idealised 
setup for the stalled shock in a supernova core and showed that 
a spherical shock is dynamically unstable to non-radial pertur- 
bations, even without neutrino heating and convection. The au- 
thors referred to this as the "standing accretion shock instabil- 
ity", or SASI, which reveals a preferr ed growth of / — 1 m ode 
defor mation and was explained by Blo ndin & Mezzac appa 
(2006) as a consequence of the propagation of sound waves 
in the volume enclosed by the shock. 

While these investigations lacked the use of a detailed de- 
scription of the neutrino ph ysics and of the equ at ion of state 
of the s upernova mediu m, [Scheck et al. (2004); Janka et al 
(200^3 1200414 12005a ). Ohrrishi et alJ feOOrfl . IBiiras et al 
(2006bJ), and Burrows et alJ ( 120061) provided results which 
demonstrate that the instability of the accretion shock also 
occurs in models which include the r elevant microphysics 
with more realism. IScheck et alJ J2004 suggested a link of 
these low-mode instabilities of the supernova shock dur- 
ing the neutrino-heatin g phase to global e xplosion asymme- 
tries (see in particular Kifonidis et al. 2006) and pulsar kicks. 
Most previous two-dimension al (2D) simulation s of success- 
ful neutrino-driven explosions ( Herant et al Jl 992llHerant et al J 
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1 1994 burrows & Frvxedl 19931 burrows et alJI 19951) failed to 
see the development of / = 1,2 modes (such an anisotropy, 
howe ver, showed up in one of the weakly exploding mod- 
els of Janka & Muilen ll996) because most of the simulations 
were done with limited computational wedges of only 90° to 
120° latitudinal width, or because very rapid explosions were 
obtained. In these cases the low-mode asymmetries were ex- 
cluded by constraining boundary conditions, or they could not 
grow in the time available between shock stagnation and re- 
vival. The latter effect may have been the reason why low-mode 
instabiliti es were not found to be do minant in the 3D simula- 
tions o flFrver & WarrerJ fcooibool . which developed explo- 
sions on rather short time scales after bounce. It is also possi- 
ble that these 3D simulations were not evolved sufficiently far 
in time to observe the formation of / = 1,2 modes. Without 
a sufficiently strong contribution of the / = 1 mode, the neu- 
tron star recoil velocities rem ain low (typically less than about 
200km/s, see lJanka & MulleJl994 . 

The main goal of the present paper (the first in a se- 
ries) is to show that global anisotropies and large NS kicks 
can be obtained naturally in the framework of the neutrino- 
driven SN explosion mechanism due to the symmetry break- 
ing by non-radial hydrodynamic instabilities, without the need 
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to resort to rapid rotation (e.g. lKotake et al, 2003), large pre- 
collapse p erturbations in the iron core jBurrows & Haves 1996; 
iGoldreic h et alJl996HLai & Goldre ichl2000l). str ong magnetic 
fields JWheeler et alJ l2002t iKotake et alJ l2004 . anisotropic 
neut rino emission associated with exotic neutr ino prop erties 
(e.g. [Fryer & Kusenko|l2006l: iFuller et al.ll2003l) . or jets dCenl 
19981 iKhokhlov et aljll999TlLai et alJl200ll) . To this end we 
present a comprehensive 2D parameter study of supernova dy- 
namics that can be considered as a significa nt improvement and 
extension of the earlier calculations of l.Tanka & Miillerl fl996) 
with respect to the treatment of neutrino transport, the assumed 
characteristics of the neutrino emission from the neutron star 
core, the inclusion of rotation, the influence of the initial seed 
perturbations, the spatial resolution, and the covered evolution- 
ary time of the supernova explosions. 

Parts of the present work were already presented in a 
Letter by IScheck et alJ ll2004l) . but a detailed description of 
both our methods and results will be given here. We proceed 
by summarising our numerical algorithms and computational 
approach in Sect. |5J and our boundary conditions and initial 
data in Sect. [5] We then give an overview of our simulations in 
Sect. 0] discussing two representative neutrino-hydrodynamic 
calculations in some detail. In Sect. [5] we explore the depen- 
dence of our simulations on the properties of the stellar progen- 
itors and on the assumed core neutrino fluxes, and establish cor- 
relations between explosion parameters and neutron star kicks. 
Section|6]is devoted to the effects of rotation. In Sect.0we re- 
turn to the neutron star recoils and investigate their robustness 
with respect to the approximations and assumptions that we 
have employed. Furthermore, we investigate the long-time evo- 
lution of the recoil velocities for a few models beyond the time 
interval of one second after core bounce, for which we have 
evolved most of our models. Estimating the terminal values of 
the NS velocities by a calibrated extrapolation procedure, we 
will speculate about the possible implications of our results for 
the velocity distribution of neutron stars in Sect. [8] A summary 
of this work and our conclusions can be found in Section|9] In 
Appendix lAl we define and tabulate some physical quantities 
of interest which characterise the different runs of our large set 
of simulations. Furthermore, we describe the post-processing 
procedures that we applied to the numerical calculations to 
compute these characteristic quantities. AppendixlBl discusses 
the solution of the hydrodynamics equations in an accelerated 
frame of reference. In Appendix ICl we analyse the explosion 
energetics of our neutrino-driven supernovae. Appendix [D] fi- 
nally details our new neutrino transport scheme. 

2. Computational approach and numerical 
methods 

2.1. Hydrodynamics and gravity 

The basic version of the com puter program t hat we employ 
for this study is described in Kifonidis et al ] d2003l) . It con- 
sists of a hydrodynamics module which is based on the direct 
Eul erian version of th e Piecewise Parabolic Method (PPM) of 
IColella & Woodward! i 19841) (augmented by the HLLE solver 
of Einfeldt 1988 to avoid the odd-even-decoupling instability), 



and a module that computes the source terms for energy and 
lepton number which enter the hydrodynamic equations due 
to neutrino absorption, scattering, and emission processes (se e 
below). The equation of state is that o fl.Tanka& Mullet! ill 996). 
In contrast to Kifonidis et al. (2003, 2006) we do not follow 
explosive nucleosynthesis in this work. This allows us to save 
a considerable amount of computer time, which is mandatory 
for carrying out an extended parameter study like the one pre- 
sented here. For the same reason, the neutrino transport in our 
simulations is described approximately (see Sect. 12. 21 and the 
dense NS core is replaced by a moving inner boundary (usu- 
ally a Lagrangian shell) whose contraction mimics the shrink- 
ing proto-neutron star. 

We include self-gravity with relativistic corrections by first 
solving the Newtonian two-dimensional P oisson equation us- 
ing a Legendre expansion according to iMuller & Steinme't3 
ill 994 . and by subsequently replacing the "spherical part" of 
the resulting gravitational potential of the 2 D mass distribu- 
tion by the "effective r elativistic potential" of iRampp & J anka 
J2002I) (for details, seefMarek et ai1 l2006l) . For describing the 
gravity of the central "point mass" (i.e., the mass enclosed 
by our inn er boundary) we use the baryonic mass where 
Eq. (53) in Ramp p & Jankal J2002I) requires the gravitational 
mass. This turned out to yield very good agreement with the 
improve d version of the ef fective relativistic potential devel- 
oped bv lMarek etaD i l2006l) . 

2.2. Neutrino transport and neutrino source terms 

The original code version of Kifonidis et al. (2003) made use 
of a simple light-bulb approximation Jjank a & Miiller 1996) 
in which luminosities of neutrinos and antineutrinos of all 
flavours were imposed at the inner boundary (which is usu- 
ally below the neutrinospheres) and kept constant with radius. 
These luminosities were typically not chosen to give accurate 
values for the fluxes prevailing below the neutrinospheric lay- 
ers, but their choice was guided by the asymptotic luminosities 
that emerge from the contracting and accreting nascent neutron 
star at large radii. This was necessary in order to cope with the 
main problem of a light-bulb approach, namely that it neglects 
changes of the neutrino fluxes and spectra that result from the 
interactions of neutrinos with the stellar matter, thus ignoring, 
for example, the contributions of the neutrino emission from 
accreted matter to the neutrino luminosity. 

In this work we considerably improve upon this former ap- 
proach by explicitly including these effects. We achieve this by 
abandoning the light bulb in favour of a gray, characteristics- 
based scheme which can approximate neutrino transport in the 
transparent and semi-transparent regimes. The approach is not 
particularly suited to handle also the regime of very large opti- 
cal depths, t. Therefore we still perform our simulations with 
an inner grid boundary at r as 10 . . . 100. However, the lu- 
minosities prescribed there have no relation to those used in 
the older light-bulb calculations. We have chosen them to re- 
produce qualitatively the evolution of the luminosities in a 
Lagrangian mass shell below the neutrinospheres as obtained 
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in recent Boltzmann transport calculations (see also Sect. 13. 1.21 
and especially Appendix lD . 2l for details). 

The transport scheme itself solves the zeroth order moment 
equation of the Boltzmann equation. The transport of neutrino 
number and energy is accounted for separately, by integrating 
two such moment equations for neutrinos and antineutrinos of 
all flavours (e,/z, t). This allows us to adopt a non-equilibrium 
description with the assumption that the spectral form is Fermi - 
Dirac, but the neutrino temperatures T Y . are not necessarily 
equal to the gas temperature T. Solving transport equations for 
neutrino number and energy, we can determine locally neutrino 
number and energy densities and thus the spectral temperatures 
T v . from the mean neutrino energies. A detailed description 
of our approximative solution of the non-equilibrium transport 
problem and the exact expressions for the employed interaction 
kernels can be found in Appendix|D] While giving qualitatively 
similar results as Boltzmann-solvers in spherical symmetry (cf. 
Sect. 14.31 . the computational cost of this approximative trans- 
port scheme is two orders of magnitude lower. 

2.3. Numerical grid and frame of reference 

We adopt 2D spherical coordinates (r, 8) and assume axisym- 
metry. Unless noted otherwise, the calculations presented in the 
following are carried out in the full sphere, i.e. for < 8 < n, 
with a grid that is equidistant in the lateral direction. A non- 
equidistant grid is employed in the radial direction whose local 
spacing, Ar, is chosen such that square-shaped cells are ob- 
tained in the convective region, i.e. Ar a rA8. Typically 400 
radial and 180 lateral zones are used. 

The outer boundary of the computational domain is typi- 
cally located at R \, ~ 2 x 10 4 km, while the inner boundary 
is placed within the forming neutron star after core bounce, at 
a Lagmngian mass shell somewhat below the electron neutri- 
nosphere. The spacing of the zones near and below the neu- 
trinospheres is chosen such that variations of the optical depth 
per zone remain smaller than a few. The baryonic matter of 
the neutron star interior to the inner boundary, M core (which is 
typically ~ 1.1 M ), is removed and its gravitational attraction 
is taken into account by assuming a point mass at r — (see 
Sect. ED- 

Since we are mainly interested in neutron star kicks in this 
paper, we need to point out that the use of the inner boundary 
condition (enclosing the NS "core") implies that the NS is at- 
tached to the centre of our computational grid. It is therefore 
not free to move relative to the ejecta during the simulation 
(unless special measures are taken, see below). This is tanta- 
mount to assuming that the NS has an infinite inertial mass. 
Two implications result from this approximation: A potential 
hydrodynamic feedback of a displacement of the NS relative to 
the ejecta is neglected, and the neutron star recoil velocity has 
to be determined indirectly in a post-processing step by mak- 
ing use of the assumption of total momentum conservation (see 
AppendixlAl. 

The relative motion between neutron star and ejecta can, 
however, be accounted for during a simulation by "wagging 
the dog", i.e. by assuming that instead of the neutron star the 



ejecta move coherently in the opposite direction of the neutron 
star's recoil. This can be achieved technically by adding the 
velocity of the relative motion to the gas velocity on the com- 
putational grid, which is tantamount to performing (after every 
time step) a Galilei transformation to a new inertial frame in 
which the neutron star core is at rest and centred at r — (see 
Appendix |B] for details). Simulations including this procedure 
will be used to investigate potential deficiencies of our standard 
assumption that the NS has an infinite inertial mass and takes 
up momentum without starting to move (see Sects. l4~4l and 171. 

3. Initial and boundary conditions 

3.1. Boundary conditions 
3.1.1. Hydrodynamics 

For solving the hydrodynamics equations reflecting boundary 
conditions are imposed at the lateral boundaries at 8 = and 
8 = 7i, while transmitting (i.e. zero gradient) boundary condi- 
tions are employed at the outer radial boundary. 

The inner boundary, which is located at the Lagrangian 
mass coordinate where we cut our initial (i.e. immediate post- 
bounce) models, is taken to be impenetrable. The contraction 
of this mass shell (and hence of the neutron star core) is mim- 
icked by moving the inner boundary of our Eulerian grid from 
its initial radius, R\ b , inwards to a final radius R f ib according to 
the expression 



Rfb(f) 



1 + (l-exp(-t/t ib ))(R\jRl-l) 



(1) 



of |janka&Mullerl Jl996). The parameter R' ib is typically in the 
range 55 km < R\ h < 85 km. 

For R l ih and f;b we use two alternative prescriptions: In 
what we henceforth will call the "standard boundary contrac- 
tion case" - as this is the or i ginal parametrization that was em- 



ployed bv I.Tanka & Mullen {1996) - we set R f ib = 15 km and 



fib = tL, where the time scale ti is connected to the luminosity 
decay and is defined in Appendix ID. 21 In the second prescrip- 
tion, the so-called "rapid boundary contraction case", we set 
R\ b = 10.5 km and r ib = 0.25 s. 

Figure ^ compares Ri\,{t) for both parameter choices with 
each other and with data fro m a supernova simulation with 
the nuclear equation of state of Latti mer & Swest vl Jl99ll) and 
with Boltzmann neutrino transport jBuras et alJl2003l) for one 
of our initial models. The standard boundary contraction results 
in a larger final radius and a slower contraction of the neu- 
tron star. The rapid boundary contraction gives results which 
are almost indistinguishable from the Boltzmann calculation. 
Although this latter parametrization is potentially more realis- 
tic than the former, we have performed the simulations with 
the "standard" case unless noted otherwise. On the one hand 
this reduces the differences compared to our previous work to 
only the treatment of the neutrino transport. On the other hand, 
the more compact neutron stars obtained with the "rapid" case 
require the use of a much finer zoning for adequately resolv- 
ing the steep density drop in the neutron star "atmosphere". 
This leads to computing times which are at least a factor of 
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Fig. 1. Evolution after core bounce of the radius corresponding 
to a mass coordinate of 1.1 M from a supernova simulation 
with Boltzmann neutrino transport (Bura s et alJl2003|) . com- 
pared to the motion of the inner boundary radius as defined 
by Eq. Q for the standard boundary contraction case (with 
£l = 1 s) and the rapid boundary contraction case. 

five larger than those of the standard case. Exclusive use of 
the rapidly contracting boundary would thus have severely re- 
duced the number of computed models, leading to much poorer 
statistics. For this reason we have chosen the rapidly contract- 
ing boundary only for a limited set of models to investigate 
whether the results obtained for the standard case change quali- 
tatively. Both cases together therefore provide information how 
the results depend on the contraction behavior of the forming 
neutron star (see Sect. 17. 3> . 

3.1.2. Neutrinos 

The boundary conditions for the neutrino emission at the in- 
ner grid boundary are chosen to be isotropic. Luminosities 
and mean energies for neutrinos and antineutrinos of all three 
flavours are imposed there in order to solve the transport prob- 
lem as described in Appendix iDl These luminosities and en- 
ergies are chosen as time-dependent functions that are con- 
strained by prescribed and varied values for the total loss of 
energy and lepton number from the core of the forming neu- 
tron star. For example, the lepton number loss during the first 
second is of order 0.1-0.2 in all our simulations, and the total 
(asymptotic) energy loss AE^^ does not exceed the gravita- 
tional energy 

E*3xl0 53 (^f (JBL-f ergs, (2) 
\M Q ) \ 10km/ g W 

which can be released during the birth of a neutron star (see 
Tables lA~TUA~5l . 

3.2. Initial models and initial perturbations 

Our calculations are started at ~ 15 - 20 ms after core bounce 
from detailed post-collapse models. We make use of four such 
models which are based on three different SN progenitors. The 



first was calculated bv lBruenr] i 19931) with a general relativis- 
tic, one-dimensional (ID), Lagrangian hydrodynamics code 
coupled to neutrino transport by multi-group, flux-limited dif- 
fusion (see h is Model WPE 1 5 LS 180). It employs the 15 M Q 
progenitor of Woos lev et a Simulations based on this 

model will henceforth be called the "B-series". 

Our second ID post-collapse model, provided by 
M. Rampp (priv. c omm.), uses a 15 M progenitor star of 
iLimongi et alJ d2000l) and was computed with the Prometheus 
PPM hydrodynamics code coupled to the Vertex multi-group 
variable Eddington factor/Boltzmann neutrino transport solver 
jRampp & Jan3l2002l) . Our "L-series" of simulations makes 
use of that model. 

We also consider two post-bou nce models that were com- 
puted for the sl5s7b2 progeni tor oflWoo slev & Weaverl (fl 995 > 
with Prometheus/Vertex by Bu ras et alJ 1120031 |20J36a) (see 
their Models sl5/lD and sl5r). The first of these models is 
from a one-dimensional simulation and gives rise to the "W- 
series" of runs, while the second is a rotating, two-dimensional 
(axisymmetric) model, which we use for our "R-series" of cal- 
culations. This latt er model has been described in detail in 
iMiilleret alJJ2004h . 

The level of numerical noise in our hydrodynamics code 
is so low that a one-dimensional, isotropic initial configura- 
tion remains isotropic, even in the presence of a convectively 
unstable stratification. Therefore we need to explicitly add 
random perturbations to trigger the growth of non-radial hy- 
drodynamic instabilities in the post-shock flow. The p ortable, 
high-quality rand om number gene rator RANLUX due to ljamesl 
i 19941 I 19961) and lLiischerl ( ll994 is employed. We apply the 
perturbation to the velocity field and typically use an amplitude 
of 0.1%. To break the equatorial s ymmetry of the rotating 2D 
model of Bur as et alJ ll2003lEo06ah . we have to add perturba- 
tions with an amplitude of several per cent, since in this model 
the initial perturbations have already grown to such a level by 
the time we map the model to our full 180° grid (see Sect. 12. 3> . 

4. Overview of the simulations 

4.1. The computed models 

Tables lA~Tl4A.5l give an overview of all our simulations (which 
were performed with the "standard boundary contraction") in 
terms of some characteristic quantities which are defined in 
AppendixlAl 

The naming convention we have chosen for the models is 
the following: The first letter denotes the initial model (i.e. 
the progenitor/post-bounce data), followed by a two-digit code 
which corresponds to the chosen value for the total asymptotic 
neutrino energy loss of the neutron star core, AE^ C0[e , in units of 
— !-q M pC 2 . Thus B18, for example, refers to a simulation based 
on the lWooslev et alJ ( fl988i)jBruenn| ( ll993l) initial data with an 
assumed release of gravitational binding energy of the core of 
AE™ cole = 0.18 M G c 2 . The second fundamental model parame- 
ter, the luminosity time scale ti, is not taken into account in the 
model names, because it has the same value for all models of 
a series. The chosen value in each case is given in the captions 
ofTables lA.lHA.5l 
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Simulations performed on a larger grid (with an outer 
boundary radius of 10 10 cm and 500 radial zones) are indicated 
by the letter "g" appended to the model name, e.g. B18-g, sim- 
ulations which account for the recoil motion of the neutron star 
contain the letter "m" in the model name, and model series 
started from different random seed perturbations are denoted 
by numbers appended to the model names. Hence Model B18- 
1 differs from Model B18 (and from Models B 18-2, B18-3 etc.) 
only in the random perturbations imposed on the initial veloc- 
ity distribution (with the perturbation amplitude being the same 
in all cases). 

Note that in Tables IA.lHA.5l the total lepton number and 
energy loss of the neutron star core, AF e core and A£ , '° < ! ore , re- 
spectively, the time-integrated energy loss in v e and v e , A.E500, 
the explosion energy, £ e xp, the anisotropy parameter, ff gas , the 
shock deformation, J s hock, the neutron star mass and recoil ve- 
locity, M ns , and v" s , respectively, as well as the correction of 
the latter due to the "neutrino rocket effect", v° s ' v , are all given 
at the time t = 1 s, at which we usually stop our simulations. 

We need to point out here that the listed neutron star ve- 
locities are not the final ones, but that even at the end of our 
simulations the neutron stars can still experience a large accel- 
eration. We therefore also give this acceleration, a° s = dv" s /df 
(averaged over the last 100 ms and without neutrino effects), 
and will attempt to estimate the final neutron star velocities in 

sect.rm 

4.2. The character of the flow 

Giving an accurate qualitative description of the flow that es- 
tablishes in our calculations is a difficult endeavour, as the 
evolution that we observe during the first ~ 300 - 400 ms is 
wildly time-dependent and extremely nonlinear. One may even 
characterise it as chaotic. The layer between the proto-neutron 
star and the supernova shock is Ledoux-unstable, because a 
negative entropy gradient is established due to neutrino heat- 
ing within ~ 50 ms after bounce. In all simulations discussed 
here, it is consequently Ledoux convection which breaks the 
initial spherical symmetry: small Rayleigh-Taylor mushrooms 
grow from the imposed random seed perturbations and start ris- 
ing towards the shock. They merge quickly and grow to fewer 
but larger bubbles that deform the shock and push it outward 
(Fig.0. 

Due to the violent motions of the rising high-entropy 
plumes the shock gets bumpy and deformed, and caustic-like 
kinks of the shock emerge where two such bubbles approach 
each other and collide. Downstream of the shock, deceler- 
ated and compressed matter forms a high-density (low-entropy) 
shell, which sits atop high-entropy material that boils vigor- 
ously as it is heated by neutrinos from below. Th e interface be- 
tween these layers is Rayleigh-Taylor unstable (Herant 1995) 
and gives therefore rise to narrow, low-entropy downflows of 
matter, which penetrate from the postshock layer to the neu- 
tron star with supersonic velocities. When they reach surround- 
ings with entropies lower than their own, the downflows are 
decelerated and their material spreads rapidly around the neu- 
tron star. The evolution of these downflows is highly dynamic. 



They form, merge with other accretion funnels, or are blown 
away by the rising buoyant matter on a time scale of 10-20 ms, 
while their number decreases with time. The most massive of 
these downflows originate from the kinks at the shock surface, 
where the deceleration of the infalling matter is weaker due to 
the (local) obliqueness of the shock. 

During this phase of violent "boiling" the shock devel- 
ops a strong, time-dependent deformation and expands slowly 
outward. In Model B12, whose evolution is shown in the se- 
quence of plots on the left side of Fig. |2] pronounced bipo- 
lar hemispheric oscillations become visible after about 150 ms. 
Such bipolar oscillations (and the consequent "sloshing") of 
the shock have been found to be typical of / = 1 mode in- 
stabilities as associated with the advective-acoustic instability 
or the SASI. Hence these instabilities are likely to dominate 
the evolution of this model in this strongly nonlinear phase. 
Note that Model B 1 2 differs from Model B 1 8 (on the right side 
of Fig. |2ji by lower neutrino luminosities at the inner bound- 
ary and, correspondingly, by a later onset of the explosion (at 
f exp = 220 ms compared to f exp = 152 ms) and a lower explo- 
sion energy of 0.37 x 10 51 erg versus 1.16 x 10 51 erg (measured 
at t — 1 s). Model B18 shows also violent convective activity, 
but no bipolar oscillations. This fact might indicate that in this 
model the convective instability might play a more important 
role. A detailed investigation of the growth of different kinds 
of non-radial instabilities in the postshock flow and their com- 
petition will be presented in a subsequent paper of this series 
(Scheck et al. 2006, in preparation). Here we only note that 
their combined effects can have a decisive impact on the explo- 
sion mechanism of supernovae, since one-dimensional coun- 
terparts of both Models B12 and B18 failed to explode. 

The highly dynamic phase of the evolution comes to an end 
around 300 - 400 ms after bounce. At that time the explosion 
is well underway, and the overall flow settles into a state of 
quasi-self-similar expansion, which is remarkab ly stable ( com- 
pare the lower panels of Fig. 13 and see also Herant 1995). Yet 
the consequences of the dynamic phase are felt long thereafter, 
since the strongly nonlinear boiling motions lead to a final mor- 
phology that is sensitive to even tiny initial differences in the 
flow. These may not only result from the influence of different 
boundary conditions, as in case of Models B12 and B18. The 
late-time morphology is even sensitive to the seed perturbation 
that we apply to trigger the non-radial instabilities. Figure [3] 
illustrates this for the B18-series of models in which the seed 
perturbation was varied as described in Sect. 14.11 It is obvious 
that the dominant mode in the flow is unpredictable. It can be 
I = 2, with two bubbles which are separated by a single accre- 
tion funnel, and which occupy roughly a hemisphere each (as 
in Model B18-4 and the original Model B18). Yet the bubbles 
may also differ significantly in size resulting in a dominance 
of the / = 1 anisotropy as in Model B18-3 (top right panel of 
Fig. |5J and Model B12 (left panels of Fig.[21>. This sensitivity to 
the seed perturbation is so extreme that the system may be de- 
scribed as exhibiting symmetry breaking in a chaotic manner. 
In fact even the same model computed on different machines 
(with slightly different 64-bit round off behaviour) may actu- 
ally end up with a different morphology. 
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Fig. 2. Entropy distributions in Model B12 (left column) and Model B 18 (right column) for different times. The figures are plotted 
such that the polar axis is oriented horizontally with "south" (6 = n) on the left and "north" (6 = 0) on the right. Dotted black 
lines mark the gain radius and white lines the supernova shock. Note that the scales differ between the plots. Convective activity 
starts with small Rayleigh-Taylor structures (f = 50 ms) which then grow and merge to larger cells and global anisotropy. In 
contrast to Model B18, the low-energy model B12 develops pronounced bipolar oscillations (compare the plots for t = 200 ms 
and t = 250 ms between both cases). After the explosion has set in, the convective pattern "freezes out" and the expansion 
continues essentially self-similarly (see the plots for t - 500 ms and t = 1000 ms). At 1 s after bounce Model B18 shows the 
emergence of an essentially spherical neutrino-driven wind expanding away from the neutron star surface (region around the 
coordinate center in the bottom right panel). 
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Fig. 3. Density distributions one second after core bounce for four simulations with the same initial and boundary conditions 
as Model B18, but different patterns of the random seed perturbations imposed on the velocity field of the initial model. The 
amplitudes of the perturbations (10~ 3 ) are the same in all cases. The morphology of the explosion depends in a chaotic way on 
the initial perturbations. 



For sufficiently high core luminosities, accretion of matter 
onto the neutron star is eventually superseded by the onset of a 
nearly spherically symmetric neutrino-driven wind (see the re- 
gion around the coordinate center in the lowermost right panel 
of Fig Eland in the left pane ls of Fig.|3] cf. also lBurrows et alJ 
1995t LTanka & Mullerl fl996). If the wind is strong enough, as 
in Model B18 where the mass-loss rate of the nascent neutron 
star by the wind is M ns = -5.1 x 10~ 2 M Q /s, it blows away the 
accretion funnel and establishes a high-entropy shell or cav- 
ity of rapidly expanding low-density material around the neu- 
tron star, which is separated from the ejecta by a strong reverse 
shock. Otherwise accretion through the funnel continues until 
more than about 1 s after bounce, as in Model B12. In this case 
the accreted material reaches infall velocities of about 1/4 of 
the speed of light, while the accretion rate at t — Is has de- 
creased to M accr « 4 x 10~ 2 M /s. Since at the same time the 
neutron star mass changes at a rate of M ns s 1.1 x 10~ 2 M Q /s, 
only a fraction of ~ 25% of the infalling matter is actually in- 
tegrated into the neutron star. The remaining 75% are heated 
and reejected with high velocity in a neutrino-driven wind that 
inflates a buoyant bubble of neutrino-heated gas in the north- 
ern hemisphere opposite to the remaining accretion funnel on 
the other side (see Fig. |2] lowermost left panel). The result- 
ing flow, which is characterised by a strong dipole mode, can 
be conveyed only incompletely with plots such as Fig. |3 and 
is much more impressively captured by movies that we have 
produced from our data 1 . 



1 A collection of movies is provided as online material of this arti- 
cle. 



These movies also show that the impact and rapid de- 
celeration of the accretion streams in the vicinity of the 
forming NS create acoustic and weak shock waves that em- 
anate fro m the neutron star's surface. As was recently also 
noted by iBurrows et alJ J2006I) . these waves propagate pre- 
dominantly into the hemisphere opposite to the accretion fun- 
nel. While traversing the low-density cavity, also the acous- 
tic waves steepen into shocks, dissipate energy, and heat the 
expanding material. However, in accordance with the results 
of LFanka & Miiller ( 1996), we observe only a modest produc- 
tion of entropy due to these waves when they propagate out- 
war d in the rapidly expanding neutrino-driven wind. In case of 
the Burrows et al. (2006) simulation the acoustic energy input 
from neutron star g-mode oscillations was found to be crucial 
for the explosion of an 11 M model. Nonspherical accretion 
was found to lead to the excitation of core g-mode oscilla- 
tions at late times (k. 300-500 ms) after bounce, whose sonic 
damping transfers a significant amount of acoustic power to 
the surrounding medium and supernova shock. G-mode oscil- 
lations are in fact also present in the outer layers of our neu- 
tron stars - i.e. in the layers that are included on our grid - 
but their amplitudes are modest and t hus they do not lea d to 
the strong consequences reported by Burrows et al .] d2006h . It 
is possible, however, that our simulations underestimate such 
effects, which would require the inclusion of the whole neutron 
star without excising the central core, and the ability to follow 
the excitation of deep modes due to a self-consistent coupling 
between accretion, core motion, and core-mode generation. On 
the other hand, our models are characterized by explosions due 
to convectively supported neutrino heating (whereas the 1 1 M 
model of Burrows et al. seemingly did not explode in that way). 
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Fig. 4. Radial profiles of the sum of the v e and v e luminosities 
for Models B12 and B18 at different times after the start of the 
simulations. 

After the explosion has been launched, our models reveal the 
development of a strong neutrino-driven wind, in which the dis- 
sipation of acoustic waves may have a different effect than in 
a more or less static gas around the oscillating neutron star. 
We point out that numerical simulations of this neutrino-driven 
outflow require very high radial resolution of the steep density 
gradient near the neutron star surface. We are not sure whether 
sufficient ly fine grid zo ning was guaranteed in the simulation 
bv lBurrows et alJ J2006). A more detailed investigation of such 
questions is in progress. 

4.3. The influence of the neutrino transport 

The fact that earlier 2D simulations, which were p erformed 
with a neutrino light- bulb description Jjanka & Mullerl [l996: 
Kifonidis et al .] l20Q3l) . were not dominated by low-order 
modes, poses the question to which extent the development of 
such global asphericity in the flow is sensitive to the treatment 
of the neutrino transport. Figure|4]shows that our new neutrino 
transport description yields radial profiles for the sum of the v e 
and v e luminosities which deviate markedly from the radius- 
independent luminosities used in a light-bulb approach: The 
luminosities are significantly modified compared to the values 
imposed at the inner boundary. After some adjustment to the 



Fig. 5. Evolution of the luminosities of v e and v e at the inner 
boundary, at the v e -sphere, and at a radius of 500 km. Note the 
different importance of the accretion contribution to the lumi- 
nosity in the low-energy explosion (Model B12) compared to 
the high-energy explosion (Model B18) and the rapid decay of 
the accretion luminosity after the onset of the explosion in the 
latter model. 



local thermodynamic conditions, which takes place in a few 
radial zones next to the inner boundary, the luminosities rise 
steeply in the cooling region below the gain radius, and decline 
slightly in the heating region farther out. The rise is caused 
by the creation of neutrinos when gravitational energy is re- 
leased during the accretion and the contraction of the neutron 
star, while the slight decline results from the absorption of the 
v e and v e in the heating region. 

The "accretion" luminosity that is produced on the grid is 
usually of the same order as the luminosity emerging from 
the core. In low-energy models, like B12, the accretion com- 
ponent exceeds the core component early on, while in high- 
energy models the core component is dominant at all times (see 
the neutrino "lightcurves" for Models B12 and B18 shown in 
Fig.HJ. 

Our transport scheme can account qualitatively well for the 
evolution of core and accretion components of the neutrino 
luminosities, for the radial and temporal evolution of the lu- 
minosities and mean energies of the radiated neutrinos, and 
for the relative sizes of the v e and v e emission (see Figs. |4j|6l 
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Fig. 6. Evolution of the mean v e and v e energies at the inner 
boundary (ib), at the v e -sphere, and at a radius of 500 km. Note 
that due to the contraction and compressional heating of the 
nascent neutron star the average energies of the radiated neu- 
trinos continue to rise until the end of our simulated evolution. 



and compare with|Lieb ejidorfer^La]j2001 iLiebendorfer et all 



2005 . iRampp & JankalbOOOl iBuras et al J [20031 l2006albl and 



Thompson et alJ l2003h . We therefore think that our current 



transport treatment is a reasonably good method for performing 
parametric explosion studies with the aim to better understand 
the role of hydrodynamic instabilities during the shock-revival 
phase of neutrino-driven supernova explosions. 

Yet, we point out here that all the previously not mod- 
elled neutrino transport effects are not the reason why the de- 
velopment of / = 1,2 modes is seen here, whereas it was 
not visible in the cal culations of I.Tanka & Miillerl ll 19961) and 
Kifoni dis et alJ d2003l) . Highly anisotropic explosion s can also 
be obtained with the li ght-bulb assumption (see Ijanka et all 
120031 1.Tanka et alj Eo04a for an example). In other words, the 
details of the functional form of L(r), which are visible in 
Fig. |2 are not decisive for the growth of the I = 1,2 modes. 
What is crucial, however, is that the explosions in the current 
models start slow ly. This was not the ca se in all but one of 
the s imulations of Ijanka & MiilleJ i 19961) and iKifonidis et all 
l l2003l) . where the neutrino luminosities were assumed to de- 
cay exponenti ally with a time scale of typically 500-700 ms 
(see Table 1 in Janka & Miiller 1996) instead of varying slowly. 



The exponential, "burst-like" decline of the neutrino light bulb 
implied fairly high initial luminosities (i.e., 4.5-5xl0 52 erg s 
for v e plus v e at the inner grid boundary) - which were required 
in case of the exponential decay for getting "typical" supernova 
explosion energies - and thus strong neutrino heating occurred 
at early times after bounce. This led to rapid explosions, which 
in turn did not leave enough time for the convective cells and 
bubbles to merge before the expansion became so fast that it 
continued in a quasi self-similar way. Since these bubbles are 
initially small, their early "freezing out" in the rapidly expand- 
ing flow had the effect that small structures (i.e. high-order 
modes) prevailed until very late times. The rapid explosions 
also caused a quick end of accretion onto the proto-neutron 
star, and therefore the neutron stars remained small. In contrast, 
the present transport description gives neutrino luminosities be- 
tween the neutrino spheres and a radius of 500 km that vary 
much less steeply than exponentially with time (see Fig. [5}. 
This leads to explosion time scales that are sufficiently long to 
allow for the formation of low-order modes. 



4.4. Acceleration and recoil of the neutron star 

As is detailed in AppendixlAl in a 2D axisymmetric calculation 
the neutron star recoil can only have a component parallel to the 
z-axis. For its calculation only the z-momenta of the gas in the 
northern and southern hemispheres need to be considered (see 
Eq. IA.20> . If the momentum density of the ejecta, p z (r, 8), is 
mirror symmetric with respect to the equatorial plane, i.e., if 



p z (r,6) = -p z (r,7T-0) 



(3) 



holds, the sum of the z-momenta of the two hemispheres van- 
ishes 



= P N + P s 

z,gas * z. 



= 0, 



(4) 



and thus also the neutron star remains at rest (cf. Eq. IA. 10i . 
The latter situation is given e.g. for an I — 2 mode, i.e. for 
two polar bubbles of equal size separated by a single downflow 
which is located in the equatorial plane. However, in general 
the expansion of the ejecta will proceed differently in the two 
hemispheres, so that a net momentum P^gas + will result. 

If a single downflow has formed, e.g., in the southern hemi- 
sphere, the expansion of the ejecta will be hampered there. 
On the other hand it will proceed unaffected in the northern 
hemisphere, and thus |P? gas | will be larger than |P? gas |. Hence, 
P z jgas will be dominated by P?„ as (which is positive since all of 
the gas in this hemisphere is moving outwards). According to 
Eq. " 



(IA. 101 . the neutron star must then be accelerated in the neg- 
ative z-direction, i.e. towards the (southern) hemisphere which 
contains the downflow. This is the situation that ultimately es- 
tablishes in Model B12 (compare Fig. [2] and Table [0, and 
which is also illustrated in the right panel of Fig. [7] In this case 
the neutron star has attained a velocity of v" s = -389km/s 
one second after core bounce, and is still accelerated with 
af = -372km/s 2 (Fig.|E}. 

Model B18 also develops a single downflow, which, how- 
ever, is located in the northern hemisphere, rather close to the 
equator. Although this model is clearly less anisotropic than 
Model B12 (which is dominated by an I = 1 mode), the 
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Fig. 7. Graphical illustration of the momentum balance between neutron star and ejecta. The largest fraction of the ejecta mass is 
concentrated in a dense shell behind the shock (bright coloured ring). For a spherical explosion (left panel) the momenta of the 
neutron star and the ejecta are zero. If the expansion of one hemisphere lags behind the other, the gas has a net momentum in 
the direction of the faster expanding hemisphere. The neutron star is always accelerated in the opposite direction, i.e. towards the 
slower moving gas (middle and right panels). This acceleration can be mediated by the gravitational attraction of the anisotropic 
ejecta (middle panel). In case accretion flows reach down to the neutron star surface (right panel), additional (hydrodynamic) 
forces may contribute, but the gravitational force, in general, remains dominant. 





Fig. 8. Left: Evolution of the neutron star velocities for Models B12 and B18. The solid curve is the neutron star recoil velocity 
derived from total momentum conservation of gas and neutron star (see Eq. lA.lQt . The dotted curve includes corrections due to 
anisotropic neutrino emission. The dashed curve is an estimate obtained by integrating Eq. (|5} over time. Right: Evolution of the 
neutron star acceleration (solid curve), as computed from the (numerical) time derivative of the solid curve shown in the velocity 
plots on the left side. Also shown are the individual terms of Eq. (|5}, corresponding to momentum transfer due to downflows, 
outflows (e.g. in the neutrino-driven wind), anisotropic pressure distribution around the neutron star, and gravitational pull of the 
anisotropic ejecta. The sum of these terms (the long-dashed curve labelled "total") agrees well with the (solid) curve obtained 
independently from total momentum conservation applied to the hydrodynamics results. 
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Table 1. Integrated momenta of the ejecta in the northern (6 < 
7r/2) and southern hemispheres as well as their sum, P Z)gas , and 
the resulting neutron star recoil velocity, v" s , at t = 1 s. 



Model 




;.gas L s J 


/• gas [1=] 


v ns [Jon-] 


B12 


1.26 x 10 41 


-0.20 x 10 41 


1.06 x 10 41 


-389.3 


B18 


1.77 x 10 41 


-3.07 x 10 41 


-1.30 x 10 41 


515.1 



larger explosion energy and faster expansion result in a larger 
|P, gas | (Table This leads to a larger neutron star kick of 
v° s = 515 km/s at t — Is, while the acceleration at this time is 
a ns _ 290 km/ s 2 (Fig.|8j. Note that these values are comparable 
to the mean pulsar birth velocities derived from observations 
(see Sect. 17. 4> . and that they are considerably higher than those 
found in earlier simulations (e.g. Janka & Miiller 1994). This is 
mainly due to the low-order modes in our calculations, which 
result in a larger gas momentum anisotropy, ar gas (cf. Eq. IA. 131 
for a definition), compared to previous work. 

The neutron star velocities shown in Fig. [8] (left panels), 
as well as their time-derivatives labelled with "derivative" and 
plotted with solid lines in the acceleration plots of Fig.|S](right 
panels), are calculated from the simulation data with our stan- 
dard post-processing approach by assuming total momentum 
conservation (see Appendix lAt. The use of this approach re- 
quires a justification, because, numerically, energy and mo- 
mentum might not be strictly conserved (i.e. up to machine 
accuracy) 2 . Moreover, momentum conservation can be guaran- 
teed analytically only if the gravitational potential can be writ- 
ten as the solution of a Poisson equation. This is, of course, 
the case for Newtonian gravity. Yet, f or the "general relativistic 
potential" of Ram pp & Janka J2002h that we used in the simu- 
lations, an equivale nt of the Poisson equation cannot be derived 
jMarek et alJ200d) . 

Since the neutron star kicks discussed in this work depend 
on the anisotropic distribution of the ejected gas, we do not ex- 
pect that the small general relativistic corrections or numerical 
errors of the mentioned kind can seriously affect the results of 
our calculations to an extent that unrealistically large values for 
the neutron star recoil velocities are obtained. This expectation 
is supported by the fact that we find similarly large neutron star 
kicks in simulations with Newtonian gravity (see Sect. 17. 3> . 
Since the NS recoil estimated from our simulations is a conse- 
quence of the anisotropic ejection of mass in the explosion, it is 
also unlikely to be linked to nonconservation of energy and/or 
momentum on a small level. In order to provide additional and 
independent evidence that the neutron star velocities estimated 
on grounds of the assumption of total momentum conservation 
are reliable, we check them by verifying the estimated neutron 
star acceleration as a sum of the different forces which con- 
tribute to a momentum transfer to the neutron star. 

For this purpose we consider a sphere of radius r$ ss 1.1 R ns 
that encloses the neutron star. The time-derivative of the neu- 
tron star momentum (and hence the neutron star acceleration 
at a certain time) can then be obtained by integrating the Euler 

2 The energy and momentum conservation properties of neutrino- 
hydrodynamic codes like the employed one are discussed in much de- 
tail in lRampp & Jankal <2002h and lMarek et alj (2006). 



equation over the volume of that sphere, resulting in 

Pm ~ ~ £ PdS-S pv v r dS + f GM ns 4 dm. (5) 

Here the individual terms account for the varying pressure f 
around the sphere, the flux of momentum through the sur- 
face of the sphere, and the gravitational acceleration due to 
the anisotropic matter distribution outside the sphere. For the 
gravitational term we assume that the matter distribution inside 
the sphere is spherically symmetric and that the gravitational 
potential is Newtonian. 

The time evolution of the acceleration corresponding to 
these terms, calculated from the data of Models B12 and B18, 
is shown in the right panels of Fig. [3] Here the second term 
has been split into two components associated with momen- 
tum flux into ("downflows") and out of the sphere ("outflows"). 
Also displayed is the sum of all terms (labelled by "total"). 
Integration over time of the latter quantity yields the dashed 
velocity curve for y" s ' syn in the left panels of Fig.[S] This should 
be compared to the solid curve (v" s ) which was computed with 
our standard post-processing approach of the gas momentum 
(and which includes the effects due to general relativistic cor- 
rections). It is evident that there are only small differences be- 
tween both results, which are significantly smaller than 10%. 
This demonstrates that the flow morphology indeed produces 
an anisotropic momentum transfer to the nascent NS, which is 
the cause of the estimated NS velocities. 

An interesting implication of Fig. [8] is the fact that the 
largest contribution to the acceleration is, in general, due to the 
gravitational term. In certain evolutionary phases also the other 
terms may contribute significantly. Yet, the total acceleration 
points nearly always in the same direction as the gravitational 
pull. Momentum transfer by pressure and gas flow (the first and 
the second term in Eq.|5} are only important as long as the inho- 
mogeneous ejecta have sonic contact with the neutron star and 
thus can exert hydrodynamic forces on the central object. This 
is the situation found in Model B18 for times before t « 0.5 s. 
After that time the supersonic neutrino-driven wind, which is 
very strong in this energetic model (due to the high neutrino lu- 
minosities) has blown away the accretion downflows from the 
neutron star. Ongoing acceleration is then exclusively caused 
by the gravitational pull of the anisotropic ejecta and decreases 
slowly as the nearly spherically symmetric wind clears the sur- 
roundings of the neutron star. Hydrodynamic forces therefore 
do not contribute at later times in Model B18. On the other 
hand, they are important at all times in Model B12. The ac- 
celeration due to the momentum flux associated with the nar- 
row downflows that reach the neutron star is usually the sec- 
ond most important term, and is directed opposite to the grav- 
itational acceleration. Anisotropics in the pressure distribution 
and wind outflow contribute on a smaller level. 

Finally, we show in Fig.[8](left) with dotted lines the neu- 
tron star velocities corrected for the effects of anisotropic neu- 
trino emission (see Appendix lAl. These effects turn out to be 
small. For Model B12 the neutron star kick is thus reduced by 
about 10%, which is unusually large. For most of our mod- 
els (including Model B18) the corrections due to anisotropic 
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Fig. 9. Dependence of some global quantities on the inner boundary luminosity. The quantities in the left column (explosion time 
scale, explosion energy, and neutron star mass) depend only on the progenitor and the boundary conditions. The quantities in the 
right column (anisotropy, neutron star velocity, and acceleration) are strongly influenced by the initial perturbations. All time- 
dependent quantities are shown at t — Is. Crosses stand for the B-series of models, stars mark results for the L-series, triangles 
denote the W-series, and diamonds refer to the R-series of models (see Sect. l3.2l for the differences between these models). 



neutrino emission are smaller than 5% (for more details, see 
Sect. ED- 

5. Dependence on the initial model and the core 
luminosity 

In this section we discuss the variation of the quantities intro- 
duced in Appendix [X] as functions of the initial model and a 
systematic variation of the imposed core neutrino luminosity 
Lib- Tables IaTTIIA. 51 give an overview. To facilitate their inter- 
pretation, we also display the most important quantities for all 
models as a function of graphically in Fig. [9] 

The results plotted in that figure show that the neutrino- 
driven mechanism as computed in our models is able to ac- 
count for different key observational aspects of supernovae and 



neutron stars simultaneously, provided that sufficient time is 
available for low-order unstable modes to form. Typical su- 
pernova explosion energies of about 10 51 erg, typical baryonic 
neutron star masses around 1.4 M@ (actually between 1.3 and 
1 .6 M Q depending on the progenitor) and high neutron star re- 
coils (with a maximum of 800km/s in Model B 18-3 after 1 s of 
post-bounce evolution, see Table lA. it . are obtained at the same 
time. 

What is also apparent is that the quantities displayed in 
Fig-Elcan be grouped in two classes, those which show a clear 
correlation with the core luminosity, Lib, and those which do 
not. Among the former are the explosion time scale, f exp , the 
explosion energy, E exp , and the neutron star mass, M ns . For a 
given initial model these integral quantities show a systematic 
variation with the boundary luminosity with only little scatter. 
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In particular, these quantities (together with the mean shock 
radius) are only weakly affected by varying the random seed 
perturbations in the way described in Sect. l4.1l (see Tables lXTi - 
IA.5i . Among the quantities which do not correlate with Ljb are 
the ones that depend on the morphology of the explosion, i.e. 
the anisotropy parameter, ar gas , the neutron star recoil velocity, 
v" s , and the neutron star acceleration, a ns . These show a strong 
sensitivity to small differences in the flow (e.g. to the initial 
perturbations), and hence essentially stochastic behaviour. The 
large scatter of neutron star recoil velocities for Models B18-1 
to B18-6 of between ~ 80km/s and 800km/s (see Table lA~fl 
illustrates this clearly. 

A higher luminosity, L^, from the neutron star core causes 
the explosion to develop faster, to become more energetic, and 
to leave behind a neutron star with a smaller mass, because less 
material can be accreted onto the core when the explosion oc- 
curs faster. The monotonic correlation between L;b and the ex- 
plosion energy E exp shows that our chosen approach to parame- 
terise our simulations can also be interpreted as one in terms of 
explosion energy. In this sense Lib and E exp can be exchanged 
as governing parameters. Note, however, that the L;b-£ e xp rela- 
tion differs between the initial models. 

A similar behaviour is also visible in Fig. ^3 f° r 
AM ga i n (f exp ), the mass contained in the gain layer at time 
f exp , as a function of Ljb for all models. In fact, it is actu- 
ally AM ga j n (f eX p) which is responsible for the progenitor depen- 
dence of the Ljb-£exp relation visible in Fig. [5] mainly because 
the recombination of free nucleons to a particles and nuclei in 
the expanding and cooling ejecta from the gain layer yields a 
significant fraction of the final explosion energy. This energy 
contribution increases with more mass in the gain layer. The 
rest of the explosion energy is due to the power of the neutrino- 
driven wind of the proto-neutron star (see AppendixO. Since 
AM ga i n (f exp ) depends on the mass accretion rate through the 
shock, there is a dependence on the density profile of the pro- 
genitor star. The different initial models reveal significant dif- 
ferences in this respect. In particular, the ILimongi et alJ pro- 
genitor exhibits considerably higher densities at the edge of the 
iron core and in the silicon shell than the Woosley et al. models, 
but this progenitor explodes later and thus at a time when the 
mass accretion rate has already decreased significantly. 

It should be noted that rotation will also affect AM ga i n (f exp ) 
(see Sect. |6j. The systematically larger mass of the gain layer 
(Fig. HUi . and the up to ~ 50% higher explosion energies of 
the rotating models compared to the non-rotating models of the 
sl5s7b2 progenitor (Fig. |9j, though, are strongly affected by 
the larger initial perturbations that we have used in the rotating 
case (see Sects. I3~2l andl6li. 

A progenitor dependence is also visible in case of f exp and 
M ns as a function of Lib, as displayed in the left column of 
Fig-El The simulations th at are based on the newer 15 M & pro- 
genitor model sl5s7b2 of IWooslev & Weaver! U995) give ex- 
plosion time scales that are systematically higher by ~ 30%, 
and final neutron s tar masses that are hig her by ~ 10% than 
those of the older Wo oslev et alJ (119881) core. On the other 
hand, the results belonging to the iLimongi et alJ J2000h pro- 
genitor again exhibit larger systematic deviations from those 
for the Woosley et al. stars. The higher mass accretion rate in 
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Fig. 10. Mass of the gain layer at the onset of the explosion 
(f exp ) as a function of the boundary luminosity for the set of 
models displayed in Fig. [5] For every initial model there exists 
an approximately linear relation between AM ga j n and L^. 



simulations with the ILimongi et alJ progenitor delays the de- 
velopment of convective motions, and thus the onset of the ex- 
plosion (fexp) compared to the other models. This prolongs the 
time that the revived bounce-shock needs to reach a certain ra- 
dius. It also reduces the explosion energy, and leads to a larger 
neutron star mass, for a given value of the boundary luminosity 
Lib- 

We focus now on the right column of Fig. [9] Recalling the 
highly nonlinear, chaotic hydrodynamic evolution in response 
to a variation of the initial perturbations described in Sect. 14.21 
one can understand that there is no clear correlation between 
Lib and the quantities a gas , v" s , and a" s , which depend on the ex- 
plosion morphology. When, however, a gas is plotted as a func- 
tion of the explosion energy (see Fig. II li . it becomes apparent 
that the area near the upper right corner in the a gas -E exp dia- 
gram, satisfying 



, / ao + Lexp / L exPi o > 1 



(6) 



with Lexp.o ~ 2 x 10 51 erg and ao ~ 0.3, is almost void. This 
indicates that high-energy explosions with large anisotropies 
are disfavoured, which is plausible because there is not suffi- 
cient time available for high-order modes to merge. In order 
to assess the impact of this result on the neutron star recoil 
by virtue of Eq. ( IA.14t . we need to consider also the scalar 

Y2\ . Figure fTTI shows 



1ATT4I 

quantity P e j, which is defined in Eq 



that it is linearly increasing with the explosion energy. Since 



OC Q^ gas -P e j 



this increase of P e j with £ exp will tend to com- 



pensate the smaller values of a gas for higher explosion ener- 
gies. Therefore high neutron star velocities (up to 800km/s at 
t — 1 s) can result for a wide range of explosion energies, or, 
equivalently, boundary luminosities (cf. Fig.|9|l. Indeed we see 
that neither bipolar oscillations nor the dominance of an / = 1 
mode are excluded when the explosion energy is moderately 
large. We expect, however, that for sufficiently large bound- 
ary luminosities the explosion time scale, and correspondingly 
Q 'gas, will become so small that the neutron star velocities will 
remain low for (very) large explosion energies. 
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Fig. 11. Anisotropy parameter a gas (upper panel) and (scalar) 
quasi-momentum of the ejecta, f e j, (lower panel, see Ea. IA.12f 
for a time of 1 second after core bounce as a function of the ex- 
plosion energy. The different symbols have the same meaning 
as in Fig. [5] 



An important result of the present work is that neutron stars 
which have attained high velocities at t — Is typically expe- 
rience very high accelerations, too (reaching up to more than 
700 km/s 2 ). This becomes apparent in the panels of Fig. 1121 
which display the acceleration at the end of our simulations 
(top) or averaged over the last half of a second, respectively, 
as a function of the neutron star recoil velocity. In fact two 
populations of models may be discriminated, a low-velocity, 
low-acceleration component and a second component extend- 
ing to much higher accelerations and velocities. The latter con- 
tains simulations with a strong contribution of the I = 1 mode, 
whereas the former is made up of models in which / = 2 or 
higher modes are dominant. Since in many of the simulations 
the accelerations are still high at t — Is, one can expect that 
their neutron star recoil velocities will significantly increase at 
still later times. We will discuss this in Sect. 17.41 



6. The effects of rotation 

We have shown that the magnitude of the neutron star recoil 
depends sensitively on the convective mode. Here we will con- 
sider the influence of rotation, which can have an effect on the 



Fig. 12. Top: Neutron star acceleration as a function of the neu- 
tron star velocity after one second. Bottom: Acceleration com- 
puted as time-averaged value over the last half of a second of 
the simulations versus neutron star velocity. The acceleration 
is multiplied by a factor <x = sign(v" s ), i.e. cr(a" s ) < corre- 
sponds to a deceleration of the neutron star. The different sym- 
bols have the same meaning as in Fig. [5] Typically, low values 
of the acceleration (<x(a" s ) < 250 km/s 2 ) are associated with 
low velocities (|v" s | < 200 km/s), while much higher values 
of cr(a" s ) are reached for higher velocities |v" s |. This suggests 
two components of the distribution, one with low velocities 
and lower average acceleration values and one with both val- 
ues being higher. The thin solid line indicates the mean values 
of <x(a" s ), binned in velocity intervals of 100 km/s. 

pattern of convection via the H0iland condition, which states 
that the flow is stable to convection if 



Ch := Cs + Ch 



(7) 



x 3 dx 



1 

+ -a 

P 



dS 



dY. 



VK \ > 0, 



IP,Y e \ u± elP,S ) 

holds (see e.g. lTassouilll978h . Here a is the total (gravitational 
and centrifugal) acceleration, and j, is the specific angular mo- 
mentum (j z = x ■ v$, where x = r sin 9 is the distance from 
the axis of rotation). In the non-rotating case the condition 
of Eq. (0 reduces to the familiar Ledoux criterion for stabil- 
ity, Cl > 0, whereas for negligible entropy- and F e -gradients 
Eq. @ becomes the Solberg-condition Cs > 0. 
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Fig. 13. Radial profiles of the specific angular momentum, j z , 
at the equator of Model R18-C for several times after the start 
of the simulation. 
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Fig. 14. Distribution of the specific angular momentum j z of 
the rotating model R18-C at t — 150 ms. Matter with larger and 
larger specific angular momentum has fallen through the shock 
(outer solid line), which leads to an overall positive gradient 
d j z 2 /dx in the gain layer. However, due to convection, which is 
suppressed only near the poles, the j, stratification and its gra- 
dient are locally perturbed. The rotation axis is oriented hori- 
zontally. 




In order to investigate how rotation changes the morphol- 
ogy, the energetics of the explosion, and the neutron star recoil 
velocities, we have computed the R-series of our models. These 
models start from a post-bounce configuration with a pertur- 
bation amplitude of several percent (cf. Sect. 13.21 . which is 
more than an order of magnitude larger than the standard per- 
turbations that we employed in our non-rotating models. Such 
a large increase of the perturbation amplitude leads to notice- 
able changes in the explosion time scale and energy. A clean 
discussion of rotationally induced effects therefore requires re- 
computing some of the non-rotating models with a higher am- 
plitude of the initial random perturbations. We do this in case 
of Models W12-C and W18-C (see Tabled), in which the same 
initial perturbations are applied as in Models R12-C and R18-C, 
whose results are also listed in Table|2] 

6.1. Evolution of the rotation rate 

The in itial rotation profile that we employ was discussed in de- 
tail by Miiller et al. (2004) (see also Fig. 1 there). Our choice 
of this angular velocity profile on the one hand maximises rota- 
tional effects in view of the most recent evolution calculations 
for magnetised rotating massive stars: it yields rotation rates 
that are more than a factor of two higher in the iron core, and 
on average a facto r of ten higher in t he silicon shell than in 
the calculations of He ger et al.l (J2004). On the other hand, it 
avoids sub-millisecond rotation of the newly formed neutron 
star, which would result for rotation rates that are significantly 
higher than the 0.5 rad/s with which our iron core was assumed 
to rotate prior to collapse. With our "standard" contraction law 
the proto-neutron star spins up due to angular momentum con- 
servation to a maximum angular velocity of about 8 x 10 2 rad/s 
at one second after core bounce. This corresponds to a rota- 
tion period of several milliseconds at the end of our simulations 
(and close to 1 -2 ms after NS contraction to a radius of 10 km). 
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Fig. 15. Entropy distribution of the rotating Model R18-C 
400 ms after the start of the simulations. The white line marks 
the supernova shock. Note the two polar downflows. The rota- 
tion axis is oriented horizontally. 

6.2. Morphology 

The rotating models evolve almost identically to the non- 
rotating ones during the first 75 ms after the start of the cal- 
culations. This is so because the Solberg-term, Cs, is negligi- 
ble in this early phase. The total angular momentum and the 
derivative of j z in the postshock region are initially rather small 
(Fig.ll3>. However, the influence of the Solberg term increases 
with time because there is a positive gradient of j z upstream of 
the shock, and matter with increasingly large specific angular 
momentum is advected into the postshock region (see Figs. 1131 
and!14>. Therefore the positive derivative of j z with x grows 
within the postshock flow. Note that, since we assume axisym- 
metry, there are no forces (other than fictitious ones) acting in <p 
direction, and hence no source terms for j z are present. The spe- 
cific angular momentum of a fluid element therefore remains 
constant, and j~ is simply carried along with the flow. 

For t > 75 ms this causes the Solberg term to become suf- 
ficiently large so that it affects the pattern of convection and 
thus leads to differences compared to the non-rotating case: All 
the rotating models develop downflows at both poles, whereas 
there is no preference for the formation of polar downflows in 
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Table 2. Rotating and non-rotating models with the same initial perturbations. For more details, see the caption of Table fA.!! 



Model 


Lib 
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fexp 


M m 




ns.v 
V z 


a T 


^gas 


^shock 
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[B] 
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[B] 


[B] 


[s] 


[M e ] 


[km/s] 


[km/s] 


[km/s 2 ] 






W12-C 


29.7 


71.5 


0.11 


68.7 


57.1 


0.40 


0.301 


1.535 


44.4 
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0.03 


0.63 


W18-C 


44.5 
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0.16 


79.0 
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0.215 
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0.21 


0.08 
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0.43 
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148.1 


0.04 


-0.03 
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Fig. 16. Radial profiles of the Solberg-term, Cs, and of the Ledoux-term, Cl, (see Eq. for = 5° ("pole") and 6 = 90° 
("equator") in Model R18-c. We show these quantities for t — 50 ms (left column) and t = 150 ms (right column). For regions 
in which Cs or Cl are negative, the absolute values are plotted as dotted lines. At t — 50 ms |OJ > ICsl and unstable regions 
(Cl + Cs < 0) are present for both latitudes. At a time of 150 ms the gradient dj z 2 /dx has become sufficiently large to make 
Cs > |ClI at the pole, and thus to stabilise the flow, whereas in the equatorial region |Csl is still small. 



the non-rotating models (compare Fig. ^] with Figs. [2] and [3j. 
These polar downflows remain stable until they are blown away 
from the vicinity of the neutron star by the neutrino-driven 
wind. The stabilisation is caused by the positive x-derivative 
of /? in the Solberg term, which is amplified by the factor 1/x 3 
near the axis of rotation. Given a positive derivative of j\, a 
matter element pushed towards the axis feels a larger centrifu- 
gal acceleration a c = jf/x 3 than the surrounding matter, and 
therefore moves back to its original position. Analogously, a 
fluid element pushed away from the axis feels a restoring force 
as well. Thus, perturbations perpendicular to the axis are sup- 



pressed and perturbations of a gas configuration in rotational 
equilibrium can only grow parallel to the axis of rotation. 

For t > 75 ms this stabilising effect of the positive angu- 
lar momentum derivative becomes sufficiently large to sup- 
press convection near the axis of rotation, i.e. to make Ch = 
Cs + Cl > there. In the rest of the postshock flow the Solberg 
term is negligible (because of its dependence on x~ 3 ) compared 
to the Ledoux term (i.e. |Csl <S |ClI) and convection is not af- 
fected much. Radial profiles of Cs and Cl illustrating this situ- 
ation are shown in Fig. 1161 

The fact that only polar downflows and no polar outflows 
form can also be easily explained. Material inside a polar 
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Fig. 17. Evolution of the shock deformation parameter d s h oc y_ 
(see Eq. IA.8> for the rotating (R-series) and the non-rotating 
models (W-series). Positive and negative values of c/ s hock char- 
acterise oblate and prolate deformation of the shock, respec- 
tively. 

downflow always consists of the lowest yVgas which is ad- 
vected through the shock (see Fig. ll4t . This guarantees a stable 
situation because the angular momentum derivative with x re- 
mains positive. In contrast, a polar outflow, i.e. a rising polar 
bubble, would contain postshock matter that would be rather 
well mixed, because a convective plume encompasses matter 
from a larger range of latitudes. Therefore such a polar bubble 
would not consist of gas with a lower j z than the infalling mate- 
rial near the poles that surrounds such a bubble. This situation 
would therefore be unstable due to the absence of a positive 
derivative d j, 2 /dx. 

Besides the differences in the pattern of convection another 
morphological difference becomes evident: The rotating mod- 
els remain more spherical, whereas the non-rotating models in 
general develop a clear prolate deformation (Fig. I17> . This is 
partly due to the polar downflows, which damp the shock ex- 
pansion near the poles. A second reason is the centrifugal ac- 
celeration of the matter between neutron star and shock. Owing 
to the accumulation of angular momentum behind the shock, 
the initially weak centrifugal forces increase, and their radial 
components reach up to 20% of the gravitational acceleration. 
Consequently the shock is pushed out farther in the equato- 
rial region than in the non-rotating models. This has interesting 
consequences for the explosion energy. 

6.3. Energetics 

In both rotating Models R12-C and R18-C that are listed in 
Table |2 the explosion energies are higher and the neutron star 
masses are correspondingly lower than in their non-rotating 
counterparts, W12-C and W18-C, also listed in that Table. 
In case of models R18-C and W18-C the energy difference 
amounts to ~ 20% (i.e. 0.2 x 10 51 erg) and must be caused 
by rotational effects. This difference builds up when the ex- 
panding and cooling neutrino-heated matter in the gain layer 
recombines from free nucleons to alpha particles (and partly 
to nuclei) and remains approximately constant in the subse- 



quent phase, in which the explosion energy increases further 
due to the neutrino-driven wind (see Fig. lC.6l and AppendixO. 
It is caused by the larger equatorial shock radius in the rotating 
model R18-C and the thus wider gain layer, which increases the 
recombining mass by 0.013 M compared to the non-rotating 
case. 

6.4. Neutron star recoil 

What are the implications of the morphological differences be- 
tween rotating and non-rotating models for the neutron star 
kicks? In the non-rotating case the highest recoil was obtained 
for Model B18-3, in which a pronounced I = 1 mode with a 
single polar downflow is present. In the rotating case such a 
flow pattern cannot establish, since we always obtain down- 
flows at both poles. However, significant asymmetries can still 
develop, since one of the polar accretion funnels may be much 
stronger than the other, or a third downflow may be dominating 
the mass distribution. High neutron star recoils are thus not pre- 
cluded, although we expect the mean and the maximum kicks 
to be somewhat smaller than in the non-rotating case. 

The results of our rather few simulations, which comprise 
only nine rotating models (see Tables lA~4l andl2l. are in agree- 
ment with this expectation: The largest neutron star recoil ve- 
locity obtained in the R-series of models is 321 km/s, whereas 
it is 640 km/s in case of the W-series (see Model W18-C in 
Table |5Jl. The average kick velocities for the R- and W-type 
models are 151 km/s, and 280 km/s, respectively. If one omits 
Model W18-C, the only W-type model with a "pure I = 1 
mode", the average kick velocity of the non-rotating models 
decreases to only 228 km/s, i.e. it is 50% larger than that of the 
rotating models. This is a relatively moderate effect if one re- 
calls that the initial angular velocity assumed in the progenitor 
core of our calculations is clearly extreme compared to the rota- 
tion rates obtained from the latest stellar evolution calculations 
llHegeretall2004h . 

6.5. Spin-kick alignment? 

Does rotation lead to an alignment of the kick direction with 
the rotation axis (the so called "spin-kick alignment")? This 
question cannot be conclusively answered on the basis of two- 
dimensional axisymmetric simulations, because in this case the 
neutron star kick is always along the rotation axis due to the 
assumed symmetry of the calculations. 

However in the context of our kick scenario also in the 
three-dimensional case effects can be imagined which may lead 
to a spin-kick alignment. On the one hand, the rotation axis is a 
preferred physical direction of the system such that the devel- 
opment of global anisotropics (e.g. polar accretion and outflow, 
bipolar oscillations) might be favoured along this direction. On 
the other hand, if the rotation period is smaller than the duration 
of the neutron star kick by a one-sided, non-axial acceleration 
(in the corotating frame), then any asymmetry will retain only 
its component parallel to the spin axis, while the perpendicular 
component will be reduced or extinguished by rotational aver- 
aging. 
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Our results seem to suggest that the first effect may be the 
more important one. For the angular momentum present in our 
models there is a tendency of anisotropies (e.g. of downflows) 
to develop preferably aligned with the rotation axis. While for 
the relatively "fast" rotation of our models (in the sense dis- 
cussed in Sects. and 16. it . the second effect may also con- 
tribute to produce spin-kick alignment, the influence of rota- 
tional averaging will be weaker for slower and possibly more 
realistic rotation. In case of "slow" rotation, i.e. for spin pe- 
riods of tens of milliseconds in the nascent neutron star and 
many hundreds of milliseconds in the neutrino-heated convec- 
tive postshock layer (which are ten times or more larger than 
in our models), rotation will be unable to enforce perfect align- 
ment of the directions of kick and spin. 

Depending on the amount of angular momentum in the 
supernova core, the hydrodynamic kick mechanism discussed 
in this paper therefore allows for both possibilities, spin-kick 
alignment for rapid neutron star rotation (f ns <k 1 s) and mis- 
alignment or incomplete alignment for long rotation periods 
(fns <; some 100 ms). This seems to be compatible with re- 
cent studies of observational constraints on neutron star kicks 
for isolated pulsa rs and for neutron stars in binary systems 
JWang et al.l20 06). although the interpretation of observations 
is still ambiguous (Johnston et al .120051) . 

7. Robustness and long-time evolution of the 
neutron star recoils 

We have seen above that rotation, even if it is noticeably faster 
than in the most recent stellar evolution models, does not pre- 
clude neutron star kicks of several hundred km/s. However, we 
have made a number of approximations in our post-processing 
analysis and used assumptions in our simulations whose impact 
on the neutron star recoil still needs to be assessed. In addition, 
we have stopped most of our simulations at a time of one sec- 
ond after core bounce, when the neutron star acceleration was, 
in many cases, still high. Hence we need to comment also on 
the later evolution of the kicks. These issues are discussed in 
the following. 

7.1. Anisotropic neutrino emission 

The neutron star recoil velocities, v" s , that are listed in 
Tables IA.lHA.5l are computed from Eq. JA.10> . i.e. they do 
not include the effects of anisotropic neutrino emission. As 
we show in Appendix|A] anisotropic neutrino emission results 
in a correction, v" s,v , of the neutron star velocity which is de- 
scribed by Eqs. (IA. 17l > and ( IA. 19i . In Sect. I4.4I we have al- 
ready seen that this correction is small for Models B12 and 
B18. This actually holds for most models. Only in a few cases 
is |v" s ' v /v° s | > 10%, and in most of these cases the neutron stars 
have small recoil velocities (cf. Tables lA~Tl4A.5l l. The correc- 
tion due to anisotropic neutrino emission in general reduces the 
kick. This can be understood from the fact that in most models 
a single prominent accretion funnel is present. The neutron star 
recoil caused by gas anisotropies is always directed towards 
this downflow, while the neutrino emission associated with the 
"hot spot" created by the downflow on the neutron star surface 



results in a "neutrino-rocket engine" that kicks the neutron star 
in the opposite dir ection (this c ircumstance was observed and 
discussed before by Frver 2004). However, the acceleration due 
to the neutrino emission remains small because the anisotropy 
parameter of the accretion luminosity is typically only a few 
per cent. 

There are several reasons for that. On the one hand, the 
neutrino-radiating tip of the accretion downstream is highly 
unstable and its position varies with time, reducing the neu- 
trino emission anisotropy by temporal averaging. On the other 
hand there is a projection factor of cos 8 of the downflow im- 
pact polar angle, 0, to be included due to the axial symmetry 
of our models. This factor also reduces the kick. Finally, the 
time scale of neutrino energy release from the accreted mat- 
ter is typically significantly longer (the cooling time scale is 
of order 10 ms) than the time scale that the gas remains com- 
pressed in the downflow tips (between 0. 1 and 1 ms) before 
it spreads around the neutron star surface. Only very close to 
the lower end of the downdrafts the density of the gas is so 
high (p i= 10 11 g cm 4 ) that the neutrino emission is extremely 
large. During its violent impact on the NS surface, the gas, 
however, overshoots equilibrium conditions. Once decelerated, 
it bounces back, reexpands immediately, and wraps around the 
neutron star at radii considerably larger than the minimal radius 
of impact. This is mainly due to the fact that the gas comes 
from far out in the progenitor star and is shock heated dur- 
ing accretion. As a consequence, its entropy is still consider- 
ably higher than the entropy of the layers around and inside 
the neutrinosphere (remember that neutrino cooling during the 
infall is too slow to cool the gas efficiently). Therefore the gas 
floats and forms an essentially spherical, high-entropy and low- 
density (p ~ 10 ll) g crrT 3 ) envelope that radiates neutrinos with 
significantly lower rates than the dense tips of the impinging 
downflows. A part of the gas is integrated in the cooling layer 
and in response to the neutrino losses settles rather slowly on 
the NS, while the other, higher-entropy part is added to the re- 
gion outside of the gain radius and is neutrino-heated until it 
is blown away again in the neutrino-driven wind. As a result, 
our models reveal that only at most 10-15% of the binding en- 
ergy of the infalling gas in the downflows are radiated highly 
anisotropically. A much larger part of the released gravitational 
binding energy is not emitted in the downflows but from the 
essentially spherical layer enwrapping the nascent NS and set- 
tling on it 3 . Due to the mass ejection in the wind, the total rate 
of energy loss in neutrinos is actually significantly smaller than 
the rate of release of gravitational binding energy correspond- 
ing to stationary accretion with the mass infall rate through the 
downflow. 



3 It should be noted that our transport approximation, which as- 
sumes that the transport equations in radial direction can be solved in- 
dependently in all angular zones of the grid, has the tendency to over- 
estimate the neutrino emission anisotropy compared to a fully multidi- 
mensional treatment. Therefore our "neutrino recoil" is likely to be an 
upper limit of the corresponding effect rather than an underestimation. 
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7.2. Inertial mass of the neutron star 

In most of our simulations we make the simplifying assump- 
tion that the inertial mass of the neutron star is infinite, i.e. the 
consequences of the neutron star motion are ignored during the 
hydrodynamic simulation. This assumption is dropped in one 
set of models which is listed in Table lA.5l In these simulations 
the feedback effect of the neutron star motion is taken into ac- 
count by changing the frame of reference in every time step, 
thus allowing the ejecta to move relative to the neutron star in- 
stead of following the neutron star motion through the ambient 
gas (see Sect. l2.3l and AppendiceslAlandlBl. 

Comparing the results obtained from both approaches for a 
sample of about 30 simulations (which made use of the bound- 
ary parameters of Models B12 and B18), one sees that any 
given model, all else being equal, develops different explosion 
asymmetry and therefore NS kick, although the explosion en- 
ergy and time scale are very similar (see Tables lA~T1 and lA3t . 
The ensemble distribution of kick velocities, however, shows 
little change, and in particular neutron star velocities in excess 
of 400km/s after 1 s of post-bounce evolution are found re- 
gardless of whether the relative motion of the neutron star is 
included or not. 

Inspecting our simulations with and without NS motion, 
we can actually not discover any obvious differences caused 
by the moving NS (the reader is invited to have a look at the 
movies for Models B12 and B12-m6 which are provided as on- 
line material of this article). We think that there are a variety 
of reasons for that. In the first place, the neutron star accelera- 
tion and velocity are typically rather small, in particular before 
and just after the explosion is launched when the acceleration 
is still unsteady (see Figs.l8landll8>. Secondly, the downflow 
deceleration and impact on the NS surface are so extremely vi- 
olent and create so much sound wave and shock activity that the 
small effect of NS motion cannot be discerned from other dy- 
namical effects. Thirdly, the downflows and also the neutrino- 
driven wind at later stages are so fast (> lOOOOkm/s) and their 
accelerations so high that the neutron star motion even with 
hundreds of km/s (but still rather modest acceleration) is only 
a small correction. 

Since the explosions in our models are triggered by neu- 
trino heating, supported by violent hydrodynamic instabilities, 
we suspect that the influence of the neutron star motion might 
just be masked and dwarfed by other dynamics so that the ex- 
plosion energy and time scale do not reveal any visible depen- 
dence. On the other hand, the nonlinear growth of the hydro- 
dynamic instabilities in the shocked layer is so chaotic that any 
small changes, independent of their detailed origin (e.g., dif- 
ferent initial seed perturbations, different rounding errors on 
different computers, different neutrino interactions, the moving 
neutron star, etc.) lead to modifications of the mass and mo- 
mentum distributions at the end of our simulations. Taking into 
account the NS motion by our transformation does not have any 
specific consequences compared to other effects that influence 
randomness. 



7.3. Neutron star contraction and gravitational 
potential 

For practical reasons, all simulations listed in Tables IA. 11 - 
IA.5I and Table |2] were performed with our "standard" pre- 
scription for the contraction of the neutron star core (see 
Sect. 13. 1 . II . although the "rapid contraction case" also dis- 
cussed in Sect. 13.1.11 is potentially more realistic. To study 
the corresponding differences, we take the "high-perturbation", 
non-rotating Model W12-C (see Sect.|5]and Table[2} as a refer- 
ence case and perform an additional simulation, Model W12F- 
c, in which we replace the slowly contracting inner boundary 
of Model W12-C with the prescription for a rapidly contracting 
proto-neutron star. Table [3] compares some quantities charac- 
terising the two models. 

Model W12F-C explodes earlier and attains a higher ex- 
plosion energy than Model W12-C. This can be explained by 
the fact that for a smaller inner boundary radius more gravi- 
tational energy is released, and that for a shorter contraction 
time scale this release occurs earlier (see also Appendix P. 
With v" s (l s) = 611 km/s the neutron star recoil velocity of 
Model W12F-C is very high. Large kicks are also found in a set 
of simulations performed with rapid boundary contraction. For 
testing this we consider for instance cases with 

1. smaller initial random velocity perturbations of 0.1% 
(Model W12F in Fig.HH}, 

2. a Newtonian gravitational potential and a constant central 
point mass chosen such that the same initial gravitational 
acceleration is o btained at a mass coordinate of 1.1 M Q as 
in the models of lBuras et all d2003l) . see Models W12F-nO, 
W12F-nl and W12F-n2 in Fig-HU 

3. a Newtonian gravitational potential and a varying central 
point mass, which is increased with time to reproduce the 
evolution of the gravitational acce leration at a mass coor- 
dinate of 1.1 M in the models of Bu ras et alJ |2003), see 
Models W12F-nv, W12F-nvl and W12F-nv2 in Fig. HE] 

All of these models have in common that they explode more 
quickly than models with the standard boundary contraction. 
Yet, for all of these variations we obtain at least one simula- 
tion with a neutron star recoil velocity of more than 400 km/s 
at t = Is (see Fig. II 8i. This demonstrates that a faster neutron 
star contraction does not preclude high neutron star kicks and 
in particular, it shows that it is not the absolute value of the 
time scale for the onset of the explosion which matters. What 
matters is the ratio of the explosion time scale to the growth 
time scale of low-mode anisotropies by hydrodynamic insta- 
bilities like convection, the acoustic-vortex cycle or the SASI 
mechanism. With the faster shrinking of the neutron star not 
only the explosion time scale decreases, but also other impor- 
tant conditions change. In particular the advection time scale 
in the postshock layer and the sound travel time between shock 
and neutron star become shorter, because the faster NS contrac- 
tion initially leads to a smaller shock radius, too. Therefore the 
velocities ahead and behind the stalled shock are higher and the 
densities in the accretion layer are larger. Since the flow pattern 
between shock and neutron star surface reacts and adjusts on a 
hydrodynamic time scale, which is significantly shorter than 
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Fig. 18. Neutron star velocities (absolute values) as functions 
of time for Models W12F-C, W12F and several other models 
with fast neutron star contraction. In six out of eight models 
the neutron star moves faster than 300 km/s at t = 1 s. 

the contraction time scale of the neutron star, the growth of 
nonradial instabilities is accelerated in case of shorter dynam- 
ical time scales in the accretion layer. Low-mode flow there- 
fore develops faster (more details will be given in Scheck et al., 
in preparation). Thus the faster NS contraction leads to larger 
global ejecta asymmetry, in spite of faster explosions. Because 
of stronger neutrino heating by the higher accretion luminosi- 
ties the explosion time scal es are indeed similarly s hort as in 
the light-bulb studies o fljanka & Miillerl dl 994 f 19961) . The pre- 
vious "burst-like" light-bulb calculations were actually rather 
disfavorable for large global anisotropies: Due to the high ini- 
tial luminosities they produced fast explosions, and because of 
the extended NS the growth of low-mode nonradial instabili- 
ties was slow. What therefore finally matters is the ratio of ex- 
plosion time scale to low-mode growth time scale, and not the 
absolute period of time in which the explosion develops. The 
faster growth of non-radial instabilities can result in even larger 
values of the anisotropy parameter a gas for "rapid" as compared 
to "standard" models with the same explosion energy. In other 
words, the envelope in the a gas - E exp plane of Fig. ITTlapDears 
shifted towards larger values of a gas for a faster contraction of 
the proto-neutron star, and hence also the average recoil veloc- 
ity (for a specified explosion energy) increases. 



In our largest sample of models sharing the same (slowly 
contracting) boundary condition, i.e. the 18 B18-like models 
listed in Tables |A~T1 and Ia31 only three simulations develop 
neutron star recoil velocities of more than 500 km/s, and only 
seven produce neutron stars with more than 300 km/ s at 1 sec- 
ond. In contrast, in just eight simulations with rapid boundary 
contraction we obtain six models with neutron star velocities 
of more than 300 km/s and three models with neutron stars 
moving faster than 500 km/s (Fig. I18l >. Better statistics would 
require more simulations, which should also be based on the 
same initial model 4 and should make use of the same gravita- 
tional potential. 

We performed some of the simulations discussed above 
with Newtonian gravity to demonstrate that the choice of the 
effective relativistic potential in our models was not essential 
for our results. We recall that only when we use the Newtonian 
gravitational potential, momentum conservation can be ex- 
pected analytically (irrespective of numerical errors and inde- 
pendent of whether the point mass is increased with time, or 
not). The results therefore show that large neutron star recoil 
velocities are not linked to any violation of total momentum 
conservation associated with the use of the effective relativistic 
potential (see the discussion in Sect. 14. 41 . 

7.4. Long-time evolution of the neutron star kicks 

In order to investigate how the neutron star recoil velocities 
evolve beyond a time of one second after core bounce, we per- 
form six exemplary long-time simulations. For these we add 
150 radial zones to our grid and place the outer grid boundary 
at a larger radius of 10 10 cm, which allows us to simulate the 
first 3-4 s of the post-bounce evolution. In three of the simu- 
lations an infinite inertial neutron star mass is assumed, while 
in the other models the hydrodynamic feedback of the neutron 
star motion is taken into account. Four of the six models are 
just continued from models which we have computed up to a 
time of one second with our standard grid. We map the corre- 
sponding data onto the larger grid at t — 750 ms and extend the 
initial model profile from the old to the larger outer boundary 
of the new grid. 

The evolution of the neutron star velocities for all of the 
long-time simulations is displayed in Fig.[l9] The neutron star 
of Model B18-3 is accelerated to more than 1200 km/s within 
3.7 s. This demonstrates that the acceleration mechanism at 
work in our calculations has the potential to explai n even the 
highe st observed pulsar velocities (see e.g. Chat teriee et alJ 
2005). The fact that Model B18-3 is the only one in our sample 

4 The comparison between B and W models is viable, however, be- 
cause both progenitor models are quite similar. 
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that produces a neutron star with more than 1000 km/s does 
not appear problematic to us. It may be a matter of low-number 
statistics and might also change when more extreme conditions 
are realized in models, e.g. by a faster contraction of the neu- 
tron star than assumed in our standard set of models. In this re- 
spect the sample of simulations plotted in Fig.ll8llooks promis- 
ing. In quite a number of those the neutron stars have large ve- 
locities at one second and also still high accelerations (see, e.g., 
Model W12F-C in TableE}. 

After 3^1 s the neutrino-driven wind has blown away all 
downflows from the neutron star vicinity and has generated a 
nearly spherically symmetric wind bubble around it. Therefore 
the neutron star acceleration diminishes and the recoil veloci- 
ties approach their terminal values. The latter can be estimated 
by extrapolating the velocities at t — Is, applying an average 
acceleration value (a" s ), as computed for the time interval be- 
tween t = 0.5 s and 1 s, over a time period Af ex trapoi, according 
to 

= y-( f = l S ) + Af extrapol x (a° s >. (8) 

The average acceleration (a" s ) is introduced as a time-average 
which is less sensitive to short-time variations and thus allows 
for a more robust extrapolation of the velocities. The factor 
Afextrapoi = 0.35 s is "calibrated" by optimising the estimates 
in case of the models of Fig. ^] The agreement of extrapo- 
lated and computed terminal velocities is typically of the or- 
der of 10%. In the following section we use Eq. (|8j to esti- 
mate the final neutron star velocities for all models listed in 
Tables |A~TMA. 51 The basic findings of our analysis do not de- 
pend on whether we use a'? s (the acceleration values at the end 
of our simulations) or (a? s ) (the mean values in the last 0.5 s) 
for extrapolating the velocities beyond the simulated period of 
one second of evolution. 



8. Implications for the neutron star velocity 
distribution 

In Sect.[5]we pointed out that Fig.^] showing the neutron star 
velocities and accelerations at t — Is, suggests the existence 
of two groups of models. One group consists of cases with low 
velocities and on average low acceleration, and the other group 
cases with high velocities and significantly higher average ac- 
celeration. The latter models are typically characterised by a 
strong / = 1 mode in the flow pattern at the end of our simula- 
tions. 

Provided the acceleration shows a trend of increasing more 
steeply than linearly with the neutron star velocity, one can ex- 
pect a growth of the separation of both populations when the 
acceleration continues over a longer period of time. Thus a bi- 
modal velocity distribution will emerge, caused by the larger 
acceleration associated with the presence of a dominant 1=1 
mode in the models of the high-velocity group. To test this pos- 
sibility, we extrapolate the neutron star motions of all of our 70 
models listed in Tables |A~TMA.5I from one second to the ex- 
pected final conditions by applying Eq. (|8j. Figurel20l displays 
both the velocity distribution at the end of the simulated evolu- 
tion (at t = 1 s; left panel) and the terminal distribution (right 
panel). 
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Fig. 19. Evolution of the neutron star velocities in six long-time 
simulations with the same boundary conditions as Model B18. 
After four seconds the acceleration has become very weak in 
all models and no significant further increase of the velocities 
is expected. For each model a thin horizontal line marks the 
extrapolated velocity value vj^ according to Eq. which is a 
rough estimate of the final neutron star velocity. 



A comparison of the panels in Fig. shows that most 
neutron stars of the high-velocity and high-acceleration group 
(which is indicated by the darker shading) accelerate to signifi- 
cantly higher velocities on time scales longer than one second. 
In contrast, only very few stars of the low- velocity group reach 
velocities in excess of 200 km/s. As a consequence, a minimum 
develops in the extrapolated distribution around 300 km/s, sep- 
arating clearly the two components in velocity space. 

We interpret this result as an interesting demonstration that 
the kick mechanism discussed here is able to produce a bimodal 
distribution of neutron star velocities simply due to the pres- 
ence or absence of a dominant / = 1 mode in the spatial distri- 
bution of the supernova ejecta. Invoking two different processes 
for neutron star acceleration is not required. It is, however, un- 
clear whether this may provide an explanation of a possible 
bimodality in the observed velocity distribution of pulsars. The 
existence of such a bimo dality is not only ambiguous, some au- 
thors finding hints (e.g. Cordes & Chernofl1ll998l iFrver et alJ 
Il998l lArzoumanian et all 120021 iBrisken et alJ 120031) . while 
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component fits differ significantly between the publications. 

Though our result is inspiring as well as tantalising, we re- 
frain from making a direct connection with observations. Such 
attempts are hampered by the limitations of our analysis, which 
does not only assume the extrapolation of Eq. (|8} to be valid 
for all cases. Our analysis is also affected by our finding that 
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Fig. 20. Histograms of the neutron star velocity distribution for the 70 models of Tables IA.lHA.5l The left panel shows the 
velocity distribution at t — Is (solid black line). The darker shaded area corresponds to the fraction of models whose neutron 
stars are moving with more than 200km/s one second after bounce. The same models are displayed with dark shading also in 
the right panel, which shows the final velocity distribution as obtained by extrapolation with Eq. 0. 



the magnitude of the neutron star kicks seems to depend on 
the neutron star contraction (see Sect. 17. 3i that is mimicked in 
our simulations by a moving inner boundary of the computa- 
tional grid. Moreover, our analysis is constrained to a set of 
15 M Q stars 5 , while linking theory with observations would re- 
quire modelling explosions for a representative distribution of 
supernova progenitors, making reasonable assumptions about 
the progenitor dependence of the explosion energy and includ- 
ing the effects from binary breakup. A large set of calculations 
would have to account for the stochastic nature of the discussed 
neutron star acceleration mechanism, thus establishing the dis- 
tribution of kick velocities as a function of the progenitor prop- 
erties. One might have the concern that in the combined data 
of all of these runs the minimum visible in the velocity distri- 
bution of Fig.[2U|is filled up. Finally, quantitatively meaningful 
calculations of neutron star kicks will ultimately have to be ob- 
tained by three-dimensional modelling. 

9. Summary and conclusions 

The aim of this work was an investigation of hydrodynamic 
instabilities in the neutrino-heated postshock layer of core- 
collapse supernovae and of the importance of such instabilities 
for the development of explosion anisotropies and neutron star 
kicks. 

For this purpose we have presented a large number of (more 
than 70) supernova simulations in two dimensions (i.e., assum- 
ing axisymmetry) for different 15M progenitor models, re- 
lying on the viability of the neutrino-driven explosion mecha- 
nism. Since this viability is still an open question and no explo- 
sions are obtained in 2D models with a detailed spectral treat- 
ment of n eutrino transport for stars more massive than about 
1 1 M (see lBuras et alJ2006albh . we triggered the explosions in 

5 The employed progenitor models, however, exhibit large differ- 
ences in core sizes and core density profiles, which actually may be 
considered as reflecting the variations over a broader range progenitor 
masses. 



our simulations by replacing the contracting core of the nascent 
neutron star by an inner boundary of the computational grid and 
assuming there suitable neutrin o luminosities from th e neutron 
star core (see the introduction in lKifonidis et al l200d for a mo- 
tivation and justification of this procedure in the light of the re- 
sults from recent Boltzmann transport supernova simulations). 
The boundary was placed at a Lagrangian mass coordinate of 
typically 1 . 1 M Q , where the neutrino optical depths were usu- 
ally 10 or higher. A systematic variation of the core neutrino 
luminosities imposed at this boundary allowed us to investigate 
the growth of hydrodynamic instabilities and the development 
of the explosion in dependence of the strength of the neutrino 
heating and thus of the size of the expl osion energy. 

In contrast to previous work Jjanka & Muileil 1 19961: 
iKifonidis et aljE003l) the neutrino luminosities of the neutron 
star core in the models presented here were not assumed to de- 
cay exponentially, but - in better agreement with transport cal- 
culations for the whole neutron star - were assumed to remain 
(roughly) constant on a Lagrangian mass shell of 1.1 M G over 
hundreds of milliseconds after bounce. With this boundary con- 
dition the approximative neutrino transport scheme developed 
for the present study ensures a radial and temporal behaviour of 
the neutrino luminosities and mean spectral energies as quali- 
tatively also found in more complete and fully consistent su- 
pernova models, i.e. the core and accretion components of the 
neutrino emission are both accounted for. 

Our main results can be summarised as follows. 

1. Random perturbations, by which we seed the growth of 
non-radial instabilities in our simulations, can grow from 
small initial amplitudes (between 0.1% and some percent 
of the fluid velocity) to global asphericities by c onvective 
instability a s well as the vortical-acoustic cycle (Foglizzo 
1200 il Eool) . provided the time until the onset of rapid 
shock expansion is sufficiently long. Once the shock ex- 
pansion gains momentum, the further growth of the insta- 
bilities, e.g. by the merging of smaller structures to larger 
ones, is quenched, and the flow pattern essentially freezes 
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out. Not the absolute time until explosion matters in this 
context, but the ratio of the explosion time scale to the typ- 
ical growth time scale of the instability. A detailed inves- 
tigation of the growth of different kinds of non-radial in- 
stabilities in the postshock flow and their competition will 
be published in a subsequent paper (Scheck et al. 2006, in 
preparation). The neutrino transport description and em- 
ployed inner boundary condition for the transport used in 
this work ensured a sufficient delay of the shock accelera- 
tion, in contrast to the ligh t-bulb p arameters employed by 
IJanka & Muilerl i 19961) and lKifonidis et alJ J2003I) . 

2. The growth of the instabilities proceeds extremely nonlin- 
early and chaotically such that the final ejecta anisotropy 
turns out to be sensitive to the initial random pattern of the 
seed perturbations as well as small differences between nu- 
merical runs (connected, e.g., to small changes in the grid 
zoning, machine roundoff errors or small differences of the 
input physics). Despite the different ejecta geometry, how- 
ever, integral parameters of the models such as the neutron 
star mass, explosion time scale or explosion energy, show 
little variability. 

3. This is different for quantities, which depend on hemi- 
spheric asymmetries. The instabilities lead to symmetry 
breaking and the ejecta can attain a net linear momentum, 
balanced by the recoil absorbed by the neutron star. In prac- 
tise, the momentum exchange was found to be mediated 
by gravitational as well as hydrodynamic forces. Typically 
the former are more important, but in cases where the neu- 
tron star accretes anisotropically over long periods of time, 
also hydrodynamic interaction can contribute significantly. 
In our standard setup for the calculations, the neutron star 
is fixed (due to the use of the inner grid boundary) at the 
centre of the grid. Since it therefore does not start moving 
in spite of momentum gain, this situation can be considered 
as a case where the neutron star is assumed to have infinite 
inertial mass. In order to test whether this affects the results, 
we performed a number of runs by imposing the negative of 
the instantaneous neutron star velocity (as calculated from 
its attained momentum) on the ambient gas on the compu- 
tational grid. This leads to a collective gas motion relative 
to the neutron star fixed at the grid centre and corresponds 
to a change of the frame of reference by applying a Galilei 
transformation after every hydrodynamics step. Of course, 
for any given model, all else being equal, the model with the 
transformation yields a different explosion asymmetry and 
a different neutron star kick (but still very similar explosion 
energy and time scale). But despite these differences of in- 
dividual simulations, we could not detect any significant 
changes of the ensemble behaviour with respect to explo- 
sion parameters or magnitude and distribution of neutron 
star kick velocities. 

4. Further tests also showed that the details of the neu- 
trino treatment, the employed gravitational potential (i.e., 
performing the simulations with Newtonia n gravity or 
an eff ective relativistic potential according to Mare k et alJ 
2006), the assumed amplitude of initial perturbations or 
the assumed contraction of the inner grid boundary (which 
mimics the shrinking of the cooling nascent neutron star) 



do not have any qualitative influence on our results for 
the neutron star kicks. Quantitatively, we discovered indi- 
cations (based on a limited set of computations, however) 
that a faster contraction of the forming neutron star - which 
may correspond to a softer nuclear equation of state or more 
rapid cooling - seems to favour higher neutron star kicks on 
average. This can be explained by a more rapid growth of 
low-mode non-radial instabilities, leading to larger values 
of the aniotropy parameter a gas for a given explosion en- 
ergy. 

5. While the neutrino flux imposed at the inner grid bound- 
ary was assumed to be isotropic in all of our simula- 
tions, the neutrino radiation at large distances from the 
neutron star could become anisotropic because of lateral 
differences in the neutrino emission and absorption. The 
biggest such differences are associated with long-lasting, 
narrow downflows through which the neutron star ac- 
cretes gas anisotropically. The gas heats up strongly upon 
falling towards the neutron star surface and getting deceler- 
ated in shocks. We found, however, that the correspond- 
ing anisotropic neutrino emission produces a neutrino- 
mediated acceleration which accounts only for small cor- 
rections to the neutron star velocities produced by the 
asymmetric mass ejection. These corrections rarely exceed 
10%. Both accelerations usually produce motion in oppo- 
site directions. The reason for this is that the neutron star 
receives a kick towards a downflow (which attracts the neu- 
tron star gravitationally or leads to a momentum deficit of 
the expanding ejecta shell on the side of the downflow), 
whereas the neutrino radiation is more intense in the hemi- 
sphere of the downflow and thus propels the neutron star in 
the other direction. Since the accretion luminosity is radi- 
ated near the neutron star surface, we do not think that our 
use of the inner boundary underestimates this effect. The 
inverse is more likely. Our radial transport tends to over- 
estimate the anisotropy of the outgoing radiation, because 
truly multi-dimensional transport would redistribute the lo- 
cally emitted neutrinos more isotropically in all directions 
instead of favourin g their ra dial propagation (see the dis- 
cussion in lLivne et all2004T) . 

6. After one second of post-bounce evolution, which was the 
period of time we simulated for most models, we obtained 
maximum neutron star velocities up to 800km/s. The mod- 
els appear grouped in two populations, one in which the 
neutron stars move with less than 200km/s and have low 
acceleration at t — Is, and another one, roughly equally 
strong, where the stars have velocities higher than 200 km/s 
and on average also higher accelerations (see Fig. [21 and 
the left panel of Fig. I20l >. The two groups differ by the ab- 
sence or presence, respectively, of a strong or dominant 
/ = 1 dipole mode in the gas distribution around the neutron 
star. The simulated models cover roughly equally a range of 
explosion energies between about 0.3 x 10 51 erg and more 
than 1 .5 x 10 51 erg. We could not detect any systematic vari- 
ations of the typical magnitude or scatter of the kick veloci- 
ties with the explosion energy. We also did not discover any 
obvious correlation of the kicks with the properties of the 
three considered 15M G progenitor stars, which exhibited 
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major differences in their core sizes and density structures. 
Rotation with a fairly high pre-collapse rate of 0.5 rad/s in 
the iron core, which in view of the most recent stellar evo- 
lution models is probably unrealistically large for ordinary 
supernovae ("see iHeger et aljE004l) . lead to slightly higher 
explosion energies (due to a larger mass in the gain layer), 
a more spherical shock surface, and the presence of down- 
flows at both poles of the rotating neutron star. This sug- 
gests a weaker contribution of an / = 1 mode in this situa- 
tion compared to the nonrotating models, where downflows 
in only one hemisphere are rather common. Although very 
extreme cases were missing in our fairly small sample of 
simulations with such rapid rotation, we nevertheless ob- 
tained kick velocities in excess of 300km/s (still rising at 
one second after bounce), and could not detect any bias to- 
wards the group with low velocities and low average accel- 
eration. 

7. The two populations in our velocity distribution at one sec- 
ond are certainly interesting in view of the possibility of 
a bimodality in the distribution of measured pulsar veloci- 
ties, which however is still controversial. We therefore at- 
tempted to derive from our set of about 70 simulations 
the distribution at the time the neutron stars have reached 
their terminal velocities. In order to do that we contin- 
ued some of our models until 3-4 seconds, at which time 
the accelerations have become very small. These models 
served for calibrating the typical period of time which a 
representative acceleration must be applied to extrapolate 
from the velocity at one second to the terminal values. 
The representative acceleration was taken as the average 
value between 0.5 s and 1 s after bounce, a choice which 
guaranteed that short-time fluctuations of the size and di- 
rection of the acceleration (which are rather frequent in 
case of low-energetic explosions) do not corrupt the ex- 
trapolation. Indeed the extrapolated velocity distribution re- 
vealed a clear bimodal structure with a minimum around 
300 km/s and a high-velocity component that extends up to 
1300km/s (right panel of Fig.l2"0l. This component consists 
of most of the neutron stars that belong to the high-velocity, 
high-acceleration group at one second. Both components 
are similar in strength, but this may depend on the choices 
of parameters for the considered set of models. The basic 
result of a bimodality, however, turned out to be very robust 
against variations of the exact way of extrapolation. 

Although the presence or absence of a pronounced I = 1 
mode in the ejecta distribution offers a natural as well as sug- 
gestive way to obtain a bimodality in the context of our hydro- 
dynamic kick mechanism, we refrain from claiming that our 
result is a strong support for the existence of such a bimodality 
in the observed distribution of pulsar velocities. There are too 
many uncertainties which might lead to a filling of the mini- 
mum of our distribution. Not only do we assume that our ex- 
trapolation law (Eq. [SJi can be applied with the same value for 
the duration of the average acceleration to all models of our 
sample, we also consider only a very constrained selection of 
progenitor models, which is not representative for the true dis- 
tribution of supernova progenitors. Though our simulations do 



not reveal systematic differences of the kicks in dependence 
of the explosion energy or progenitor structure, we do not feel 
able to exclude that a correlation of both over the range of su- 
pernova progenitors could conspire such that the bimodality of 
Fig-EHlgets wiped out. 

The proposed hydrodynamic kick mechanism, however, 
leads to an unambiguous prediction, which might be tested by 
future detailed observations of supernova remnants: The mea- 
sured neutron star velocity should be directed opposite to the 
momentum of the gaseous supernova ejecta. This is different 
from many theories which explain pulsar kicks by anisotropic 
neutrino emission from the nascent neutron star. In that case 
the direction of the acceleration can be independent of ejecta 
asymmetries. 

Apart from all the assumptions and approximations enter- 
ing this work and discussed in detail above, the biggest defi- 
ciency of the present analysis is the fact that it is based on simu- 
lations which assume axial symmetry with the polar axis being 
a coordinate singularity that is impenetrable for the fluid flow. 
Currently it is neither clear to which degree pronounced I - 1 
modes of the ejecta distribution and long-lasting downflows of 
matter to the neutron star can develop in the three-dimensional 
environment, and how common they are, although first 3D sim- 
ulations with the setup and input physics described here are 
promisin g (they will b e presented in a future publication, but 
see Janka et alJ200 5b for some results). Nor is it clear what the 
distribution of neutron star recoil velocities from 3D models 
will be. The large number of long-time simulations required by 
the stochastic nature and long duration of the proposed hydro- 
dynamic kick mechanism, is currently out of reach due to its 
enormous demand of computer time. Our results must there- 
fore be considered as indicative but they are far from providing 
definitive answers. 
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Appendix A: Post-processing of the simulations 

In the following we define and tabulate some interesting char- 
acteristic quantities that were evaluated for our about 80 hydro- 
dynamic models by post-processing the data of the simulations. 
To keep the evaluation as straightforward as possible we some- 
times employ approximations which we will detail below. 

The inner boundary condition for the neutrinos is con- 
strained by the parameters ti, AE l ° l cors Xt) and AF e ,core(0- Their 
definitions are given in Appendix lD.21 while their actual values 
(at t — 1 s in case of the time-dependent quantities) are listed in 
Tables IA.lHA.5l A characterising value for the neutrino lumi- 
nosities imposed at this inner boundary is 

L ib (f) = L e ^(R ib , t) + L e ,y c (R ib , t), (A. 1) 

where L e Ve and L e i > e are the energy luminosities of v e and i> e 
defined in Appendix lD.il Equation flA. lb neglects the contri- 
bution from heavy-lepton neutrinos, whose interactions in the 
computational domain are less important than those of v e and 
v e , and who, in particular, do not contribute to the neutrino 
heating behind the shock at a significant level. 

We also consider the sum of the v e and v e luminosities at a 
radius of 500 km, 

L 500 (t) = L e , Vo (r = 500 km, f) + L e ^{r = 500 km, f), (A.2) 

and define the time average of this quantity in the time interval 
[0, ? exp ] as 

(£ 5 00> = fexp~' I " £500(0 df. (A.3) 

Jo 

The value of (L500) represents approximatively the (v e + v e ) 
luminosity that is responsible for the energy deposition behind 
the supernova shock until the explosion sets in at a post-bounce 
time t = f exp . Therefore the difference between (L500) and 
can be considered as a rough measure for the radial change 
of the neutrino luminosities in contr ast to their constancy in 
case of the light-bulb scheme used bv ljanka & Mullerl (H996). 
Furthermore, we list the total energy in v e and v e neutrinos that 
streams through a sphere with a radius of 500 km in the time 
interval [0, f], 

A£ 5 oo(0= f CL e ,v e +£ e ,v e )(r = 500km,r')dr'. (A.4) 
Jo 

The explosion energy, E exp , of a model is defined as the sum 
of the total energy of all zones of the grid where this energy is 
positive, i.e. 

-Eexp = etot <'' Am ''' ( A ' 5 ) 

£",„,,,- >0 
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Table A.l. Simulations based on the Wooslev et alJ fl988)/Bru ennl dl993h post-bounce model. The luminosity time scale ti is 
1 s. Unless noted otherwise the inertial mass of the neutron star is assumed to be infinite for these and the simulations listed in the 
following tables, i.e. the neutron star takes up momentum but cannot move on the grid. For the definitions of the listed quantities 
see the main text. All time-dependent quantities are given at a time t = 1 s, when we terminated the simulations. Energies are 
given in units of 1 B = 1 Bethe = 10 51 erg. 



Model 






AY 

'- w e,core 


(£-500 ) 


A.E500 


^exp 


*iexp 


M m 


v" s 




a T 




^shock 




[B/s] 


[B] 




[B/s] 


[B] 


[B] 


[s] 


[Mo] 


[km/s] 


[km/s] 


[km/s 2 ~] 






BIO 


24.7 


59.6 


0.09 


57.1 


45.9 


0.19 


0.294 


1.426 


-164.1 


44.4 


-180.2 


0.24 


0.67 


Bll 


27.2 


65.5 


0.10 


58.8 


46.3 


0.27 


0.280 


1.401 


-23.6 


0.7 


-248.9 


0.03 


0.97 


B12 


29.7 


71.5 


0.11 


60.6 


48.7 


0.37 


0.220 


1.399 


-389.5 


45.0 


-372.4 


0.32 


0.06 


B12-1 


29.7 


71.5 


0.11 


60.5 


47.5 


0.33 


0.228 


1.377 


72.8 


-4.7 


47.9 


0.07 


0.22 


B12-2 


29.7 


71.5 


0.11 


60.9 


48.5 


0.39 


0.212 


1.391 


85.8 


9.7 


345.7 


0.07 


0.82 


B12-3 


29.7 


71.5 


0.11 


60.9 


46.5 


0.38 


0.207 


1.369 


242.0 


2.0 


464.3 


0.18 


0.97 


B12-4 


29.7 


71.5 


0.11 


61.1 


47.7 


0.35 


0.216 


1.385 


-115.1 


20.4 


-154.2 


0.10 


0.51 


B12-5 


29.7 


71.5 


0.11 


61.0 


47.8 


0.33 


0.211 


1.387 


-206.9 


11.6 


-483.1 


0.19 


0.52 


B13 


32.2 


77.5 


0.12 


62.4 


49.6 


0.45 


0.188 


1.378 


-355.3 


32.0 


-408.0 


0.25 


0.36 


B14 


34.6 


83.4 


0.13 


63.6 


49.6 


0.51 


0.198 


1.345 


-128.0 


-11.2 


-66.7 


0.07 


0.40 


B15 


37.1 


89.4 


0.14 


65.3 


50.3 


0.65 


0.162 


1.318 


36.1 


-1.0 


36.0 


0.02 


0.27 


B16 


39.6 


95.3 


0.15 


66.3 


51.8 


0.81 


0.160 


1.305 


-214.6 


-2.6 


-334.4 


0.08 


0.57 


B17 


42.1 


101.3 


0.15 


67.6 


53.3 


0.95 


0.146 


1.289 


-25.5 


14.8 


-102.6 


0.01 


0.05 


B17-1 


42.1 


101.3 


0.15 


67.8 


53.4 


0.92 


0.160 


1.290 


-354.0 


5.6 


-202.2 


0.12 


0.31 


B18 


44.5 


107.3 


0.16 


68.3 


54.8 


1.16 


0.152 


1.275 


515.3 


5.2 


290.5 


0.15 


0.42 


B18-1 


44.5 


107.3 


0.16 


68.4 


54.7 


1.12 


0.154 


1.274 


-126.5 


-0.8 


-49.1 


0.04 


0.20 


B18-2 


44.5 


107.3 


0.16 


68.9 


54.7 


1.14 


0.152 


1.268 


82.5 


-5.2 


16.5 


0.02 


0.07 


B18-3 


44.5 


107.3 


0.16 


68.8 


57.1 


1.15 


0.142 


1.305 


798.8 


-41.2 


552.1 


0.24 


-0.06 


B18-4 


44.5 


107.3 


0.16 


68.2 


54.6 


1.14 


0.150 


1.272 


-171.6 


4.0 


65.7 


0.05 


0.46 


B18-5 


44.5 


107.3 


0.16 


68.5 


55.2 


1.09 


0.164 


1.280 


-121.8 


-0.9 


15.4 


0.04 


-0.02 


B18-6 


44.5 


107.3 


0.16 


68.7 


55.4 


1.11 


0.160 


1.283 


502.1 


-20.6 


220.0 


0.15 


-0.06 


B18-gl 


44.5 


107.3 


0.16 


68.7 


54.5 


1.12 


0.142 


1.269 


-60.3 


3.9 


-55.4 


0.02 


0.06 


B18-g2 


44.5 


107.3 


0.16 


68.7 


54.8 


1.12 


0.138 


1.273 


267.9 


-8.1 


126.7 


0.08 


0.28 


B18-g3 


44.5 


107.3 


0.16 


68.5 


54.9 


1.10 


0.150 


1.274 


-7.4 


-3.5 


0.9 


0.00 


0.02 


B18-g4 


44.5 


107.3 


0.16 


68.7 


54.5 


1.16 


0.132 


1.270 


-416.8 


1.7 


-150.9 


0.11 


0.37 


B19-gl 


47.0 


113.2 


0.17 


69.6 


55.9 


1.31 


0.148 


1.253 


-273.8 


0.3 


-96.7 


0.07 


0.41 


B19-g2 


47.0 


113.2 


0.17 


69.5 


56.0 


1.33 


0.148 


1.255 


188.5 


6.4 


48.8 


0.05 


0.15 


B19-g3 


47.0 


113.2 


0.17 


70.0 


56.6 


1.26 


0.132 


1.263 


366.6 


1.1 


183.7 


0.10 


0.13 


B19-g4 


47.0 


113.2 


0.17 


70.0 


56.8 


1.33 


0.130 


1.267 


477.1 


-18.3 


195.6 


0.12 


-0.02 


B20 


49.5 


119.2 


0.18 


71.0 


57.3 


1.49 


0.128 


1.238 


133.2 


5.6 


52.6 


0.03 


0.40 


B21 


51.9 


125.1 


0.19 


72.1 


58.5 


1.72 


0.122 


1.222 


30.6 


-0.9 


-20.2 


0.01 


0.24 



Table A.2. Simulations based on the Limon gi et alJ JioOOt/Rampp post-bounce model. The luminosity time scale ti is 0.7 s for 
these simulations. For more details, see the caption of Table TA.ll 



Model 


L\b 




Ay 


(-^500 > 


A.E500 


F 


fexp 


M ns 


V T 


ns,v 
V- 


af 




^shock 




[B/s] 


[B] 




[B/s] 


[B] 


[B] 


[s] 


[M ] 


[km/s] 


[km/s] 


[km/s 5 ] 






L12 


42.4 


94.6 


0.13 


90.7 


70.7 


0.51 


0.321 


1.677 


278.5 


-12.9 


334.3 


0.24 


0.11 


L13 


45.9 


102.5 


0.14 


91.7 


69.2 


0.68 


0.268 


1.620 


-92.6 


-5.9 


-333.6 


0.05 


0.77 


L14 


49.5 


110.4 


0.15 


94.6 


72.8 


0.81 


0.280 


1.628 


482.1 


-22.0 


297.1 


0.26 


0.31 


L15 


53.0 


118.3 


0.17 


96.2 


75.2 


1.02 


0.266 


1.617 


-239.5 


-3.9 


-378.5 


0.10 


0.63 


L16 


56.5 


126.2 


0.18 


97.8 


76.3 


1.07 


0.256 


1.586 


-437.9 


12.8 


-715.2 


0.17 


0.47 


L17 


60.1 


134.0 


0.19 


100.3 


77.4 


1.19 


0.256 


1.558 


-24.7 


5.5 


-47.6 


0.01 


0.37 



where i is the zone counter, Am, the mass contained in zone 
2, and the total specific energy e tot is given by the sum of the 
specific gravitational, kinetic, and internal energies, 



For the sake of simplicity we use here the one-dimensional 
Newtonian expression 



e to t 



Sgrav + 2 V + e ' ml - 



(A.6) 



v(r) = - 



GM(r) 



(A.7) 
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Table A.3. Simulations based on the non-rotatin dWooslev & Weaverl d^QShlBuras et al.l (2003) post-bounce model. The lumi- 
nosity time scale ti is 1 s for these simulations. For more details, see the caption of Table TA.ll 



Model 


— r — 




Ay 


"77 — T~ 




— t. 

^exp 


'exp 






v" s ' v 


a m 


CY ...... 


'■'MlOCk 




[B/s] 


[B] 




[B/s] 


[B] 


[B] 


[s] 


[M ] 


[km/s] 


[km/s] 


[km/s 2 ] 






W10 


24.7 


59.6 


0.09 


64.3 


55.4 


0.21 


0.420 


1.568 


-129.8 


42.1 


-443.1 


0.15 


0.81 


W12 


29.7 


71.5 


0.11 


69.0 


53.9 


0.31 


0.322 


1.501 


-97.7 


-9.7 


-132.5 


0.10 


0.61 


W12-1 


29.7 


71.5 


0.11 


68.0 


59.5 


0.32 


0.374 


1.563 


-363.8 


81.2 


-377.0 


0.32 


0.13 


W14 


34.6 


83.4 


0.13 


72.9 


56.6 


0.46 


0.250 


1.473 


-62.0 


-1.5 


66.1 


0.04 


0.37 


W16 


39.6 


95.3 


0.15 


76.0 


58.5 


0.67 


0.244 


1.430 


287.2 


-5.5 


464.2 


0.14 


0.68 


W18 


44.5 


107.3 


0.16 


79.3 


61.5 


0.89 


0.226 


1.401 


-283.6 


4.2 


-290.1 


0.11 


0.44 


W20 


49.5 


119.2 


0.18 


82.0 


63.5 


1.36 


0.216 


1.354 


-377.3 


0.6 


-277.0 


0.10 


0.39 



Table A.4. Simulations based on the rotating lWooslev & Weaver! \ 1995l)lBuras et aD (120 03 ) post-bounce model. The luminosity 
time scale ti is 1 s. For more details, see the caption of Table fA.ll 



Model 


Lib 


A plot 


AY 


(L500 > 


AEsqq 


^exp 


fexp 


M m 




ns,v 
V z 


a T 


ttgas 


^shock 




[B/s] 


[B] 




[B/s] 


[B] 


[B] 


[s] 


[M e ] 


[km/s] 


[km/s] 


[km/s 2 ] 






R10 


24.7 


59.6 


0.09 


59.9 


48.8 


0.25 


0.418 


1.521 


-15.4 


-14.3 


-118.7 


0.02 


-0.02 


R12 


29.7 


71.5 


0.11 


64.6 


49.9 


0.50 


0.316 


1.461 


-235.8 


17.5 


-203.4 


0.16 


0.15 


R14 


34.6 


83.4 


0.13 


69.2 


52.4 


0.69 


0.264 


1.420 


88.4 


14.6 


86.9 


0.04 


0.15 


R16 


39.6 


95.3 


0.14 


71.9 


56.0 


0.98 


0.256 


1.396 


321.2 


-8.9 


210.1 


0.11 


0.06 


R18 


44.5 


107.3 


0.16 


75.8 


58.3 


1.24 


0.232 


1.349 


-4.8 


-3.7 


-26.7 


0.00 


-0.07 


R18-g 


44.5 


107.3 


0.16 


75.8 


58.5 


1.23 


0.226 


1.352 


-113.9 


2.1 


-188.1 


0.03 


0.07 


R20 


49.5 


119.2 


0.18 


78.8 


60.9 


1.64 


0.214 


1.309 


280.1 


0.8 


123.9 


0.06 


0.14 



Table A.5. Simulations based on the Woosl ev" et al.l fl988),Bru ennl lll993h post-bounce model. The luminosity time scale ti is 
1 s. For more details, see the caption of Table IA.fi Different from the models listed in all other tables, the recoil motion of the 
neutron star was accounted for in the simulations listed here (as described in Sect. l2.3l and AppendixIBt. 



Model 


Lib 




AV 

* e.core 


(L500) 


A..E500 


-^exp 




M m 


vf 


ns.v 
V z 


a T 




<4hock 




[B/s] 


[B] 




[B/s] 


[B] 


[B] 


[s] 


[M e ] 


[km/s] 


[km/s] 


[km/s 2 ] 






B12-ml 


29.7 


71.5 


0.11 


60.9 


47.4 


0.36 


0.226 


1.384 


-56.8 


-1.7 


-208.2 


0.06 


0.48 


B12-m2 


29.7 


71.5 


0.11 


60.9 


47.7 


0.31 


0.222 


1.385 


-100.0 


19.1 


-63.5 


0.10 


0.72 


B12-m3 


29.7 


71.5 


0.11 


61.2 


47.8 


0.38 


0.210 


1.388 


272.6 


-16.5 


91.9 


0.23 


0.35 


B12-m4 


29.7 


71.5 


0.11 


60.9 


47.0 


0.35 


0.209 


1.378 


-104.3 


-7.4 


-197.2 


0.09 


0.43 


B12-m5 


29.7 


71.5 


0.11 


60.8 


47.9 


0.35 


0.219 


1.389 


365.6 


-10.1 


219.1 


0.32 


0.47 


B12-m6 


29.7 


71.5 


0.11 


60.7 


48.4 


0.36 


0.229 


1.395 


-334.1 


42.4 


-462.9 


0.30 


0.26 


B18-ml 


44.5 


107.3 


0.16 


68.9 


54.9 


1.12 


0.136 


1.274 


43.3 


-4.8 


-108.8 


0.02 


0.12 


B18-m2 


44.5 


107.3 


0.16 


68.9 


54.8 


1.14 


0.139 


1.273 


-86.8 


-1.1 


-31.1 


0.03 


0.20 


B18-m3 


44.5 


107.3 


0.16 


68.8 


55.3 


1.12 


0.131 


1.281 


76.4 


-8.8 


-11.4 


0.03 


0.39 


B18-m4 


44.5 


107.3 


0.16 


68.5 


54.9 


1.14 


0.150 


1.274 


-118.7 


14.5 


-156.4 


0.05 


0.13 


B18-m5 


44.5 


107.3 


0.16 


68.3 


54.7 


1.12 


0.166 


1.273 


-339.7 


-4.5 


-152.4 


0.13 


-0.06 


B18-m6 


44.5 


107.3 


0.16 


68.6 


55.4 


1.12 


0.166 


1.283 


-439.3 


14.0 


-194.5 


0.17 


0.04 


B18-m7 


44.5 


107.3 


0.16 


68.8 


54.7 


1.12 


0.138 


1.272 


109.2 


8.6 


2.1 


0.04 


0.38 


B18-m8 


44.5 


107.3 


0.16 


69.3 


54.5 


1.13 


0.134 


1.269 


455.0 


-4.1 


187.4 


0.17 


0.05 



to evaluate the gravitational energy, neglecting the relatively 
small general relativistic corrections, which have been taken 
into account in the simulations. 

The explosion time scale, f exp , is defined as the time after 
the start of the simulation when E exp exceeds 10 48 erg. It turns 
out that the exact choice of this threshold value does not matter 
very much. Other definitions of the explosion time scale (e.g., 
linked to the time when the expansion velocity of the shock ex- 
ceeds a certain value) do also not lead to qualitatively different 
results. 



To characterise the deviation of the shape of the supernova 
shock from a sphere we introduce a shock deformation param- 
eter, 

_ max (RM cos 0) - min (R s (0) cos(fl)) 
shock -" 2xmaxOR s (0)sinr3) l ' (A ' 8) 

where R s (6) is the local shock radius as a function of polar an- 
gle 9. The numerator and denominator in Eq. ( IA.8> are the max- 
imum shock diameters in projection on the symmetry axis and 
perpendicular to it, respectively. A prolate deformation leads to 
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a positive value of d S hock, an oblate deformation gives a nega- 
tive value. Note that a linear shift of the shock surface in the 
direction of the z-axis does not change <i s hock- 

The neutron star mass and the neutron star radius are con- 
sidered to be associated with a certain value of the density, 
p ns = 10 n g cirT 3 . The neutron star radius, R as , is then simply 
defined as the radius where the lateral average of the density is 
equal to p ns , and the baryonic mass of the neutron star, M ns , is 
given by the sum of the central point mass and the mass integral 
over all grid zones with densities > p ns . 

In evaluating the neutron star recoil velocity, v ns , we have to 
distinguish between simulations in which we consider the neu- 
tron star motion relative to the ejecta by changing the frame of 
reference after each time step (see Sect. l2.3l and Appendix 151. 
and simulations in which this motion is not accounted for. In 
the first case no post-processing is required, because the neu- 
tron star velocity is given at all times by the accumulated effects 
of the Galilei transformations applied until time t or time step 
m, 

M0 = J] A< ore , (A.9) 
«=l,.„,m 

where Av" ore is given by Eq. ( IB.9i . In the second case, v ns is 
computed a posteriori, by making use of linear momentum con- 
servation. The total momentum of the system, i.e. the sum of 



The neutron star acceleration corresponding to the velocity 
change at a given time is calculated by finite differences: 



Jn) 



„(«+!) 



„(«-!) 



f(«+l) _ f(n-l) 



(A.15) 



In computing the recoil velocity according to Eqs. (IA.10> 
and jA.14> . we have so far neglected the fact that the neutron 
star may also be accelerated by anisotropic neutrino emission. 
While our core luminosities at the inner grid boundary are as- 
sumed to be isotropic at all times and no neutron star accel- 
eration can result from these, direction-dependent variations 
of the thermodynamic variables in layers close to the neutron 
star surface develop during the simulations and ultimately lead 
to anisotropies of the neutrinospheric emission of neutrinos. 
In particular, density inhomogeneities and local hot-spots (in 
temperature) occur as a consequence of narrow accretion flows 
that transport gas from the postshock layers to the neutron star, 
where they are decelerated in shocks and radiate away energy 
in neutrinos. The anisotropy of this neutrino emission can give 
rise to a "neutrino rocket effect", whose magnitude can be esti- 
mated by considering the integrated momentum of the escaping 
neutrinos. 

For a transport scheme along radial rays like ours, the neu- 
trino momentum density has only a radial component and can 
thus be written as (see also AppendixlDl 



of the surrounding gas on the computational grid, Pgas, is ini- 
tially zero (because all models that we consider are spherically 
symmetric or equatorially and axially symmetric just after col- 
lapse). Hence we have for all times 



MO = ~Pgas(t)/M m (t), 



(A. 10) 



and v ns (f) can be determined by evaluating the neutron star 
mass and the momentum integral of the ejecta gas, 



Pg-as(t) 



J 

JR n .<r<o 



pv dV. 



Here d V = r 2 sin dr d6 d(f>. Note that the volume integral in 
Eq. JA.1U is limited by the outer boundary of our Eulerian 
grid and that the momentum flux associated with anisotropic 
mass flow over the grid boundary would have to be taken into 
account. 



Equation JA.lOt may actually also be coined in terms of 
an an i sotropy param eter of the ejecta, a gas (see Janka & Muller 
1 1994 Herant 1995). To accomplish this, we make use of the 
following scalar quantity 



p\v\ dV, 



(A. 12) 



which has the dimension of a momentum. Then we can write 
the anisotropy parameter as 

a g zs := I-Pgasl / Pej, (A. 13) 

and the absolute value of the neutron star velocity as 

I Vns I = Cgas Ay / M DS . (A. 14) 



p v e r = 



■ 6 r — _ 

C 1 



(A.16) 



where F v is the local neutrino energy flux and e r the unit vector 
in the radial direction. The integrated neutrino momentum at 
time t is then given by 



Py(t) 



X ib < 



p v e r dV 



p v e r dV + 



f dt£ 

Jo Jr=fi ol 



p v ce r dS, (A. 17) 



(A. 11) with the surface element dS = r 2 sin# dO d<p. Here the sur- 
face integral accounts for the fact that a significant amount of 
neutrino momentum may have left our grid through the outer 
boundary by the time t. The momentum of the neutron star, in- 
cluding now also the effect of anisotropic neutrino emission, is 



P™ = -(P m + P v ), (A. 18) 

so that the neutron star velocity, corrected for the recoil by 
anisotropic neutrino emission, can be written as 

Vns.cniT — V n c + V 



nsv = -P m IM m - P v /M n 



(A. 19) 



We finally note that for symmetry reasons P gas and P v , and 
thus also P ns and v ns , can have only a component parallel to 
the symmetry axis, i.e. along the z-axis, in 2D axisymmetric 
calculations. Equation ( IA.1U . for instance, therefore reduces 
to 



z.gas 



In 



2n 



dr 

R n , JO 

Jr»oo p. 
dr 



= P* 



d0r z sin6» p z (r,0) 

dd r 2 sin [ p,(r, 0) + p z (r, n - 0) ] 

(A.20) 
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Here p z (r, 9) = p (v r cos 8 - Vg sin 0) is the z-component of the 
momentum density of the gas, and ^ gas and P? gas are intro- 
duced as the z-momenta of the gas in the northern and southern 
hemispheres, respectively. 



Appendix B: Hydrodynamics in an accelerated 
frame of reference 

In an inertial frame of reference the hydrodynamic equations 
are given by 



dp 

! + V.( P v) = 0, 



p|j£ + („. V)v] + W> = P g. 

-((pE + P) v) = vpg, 

ot 



(B.l) 
(B.2) 
(B.3) 



where p is the density, v is the velocity, P is the pressure, g is 
the gravitational acceleration and E - e + v 2 /2 is the sum of 
internal energy, e, and kinetic energy, e^ n , per unit mass. 

Let AF be a frame of reference that coincides with an iner- 
tial frame IF at time t — and accelerates with a constant rate 
a in z-direction, a = ae z . The Cartesian coordinates of both 
frames are then related by 



(x',y',z, t') - (x,y,z- at 1 12, t) 



(B.4) 



(primed quantities are used for the accelerated frame), which 
implies that 

dz(x,y,z,t)/dt = -at and dz{x ,y' ,z ,t')/dt' = at. (B.5) 

For density, pressure, velocity, kinetic energy and gravitational 
acceleration the following relations hold: 



Similarly, it can be shown that the additional terms arising 
in the energy equation JB.3I for an accelerated frame of ref- 
erence are of order f 2 and can also be neglected, as long as 
\at\ <k |v.| holds. 

Within a typical time step At of a supernova simulation (of 
order 1CT 6 s) the condition \aAt\ «: |v z | is satisfied, because the 
maximum neutron star accelerations are of 6>(10 8 cm/s 2 ), and 
hence |aAf| = (9(100 cm/s), which is much smaller than the 
relevant velocities in the simulations, which are of (9(10 6 cm/s). 
Thus a solution of the inertial frame hydrodynamics equations 
with the simple replacement g — > g' — g - a should yield an 
excellent approximation to the solution of the hydrodynamic 
equations in the accelerated frame. 

Unfortunately, in the present problem the neutron star ac- 
celeration, and hence the instantaneous acceleration of the 
frame, a(t), is not known a priori, because it is coupled to 
the solution of the hydrodynamic problem during a consid- 
ered time step. Therefore we need to make use of an operator- 
splitting approach, in which we first ignore the acceleration of 
the frame of reference and simply solve the inertial frame hy- 
drodynamics equations (just using the gravitational accelera- 
tion g). We can then compute the current value of a(t), which 
is assumed to be constant over the time step, using momentum 
conservation: The sum of the momenta of the neutron star core, 
/"core, and the matter on the numerical grid, Pgnd> is conserved 
and initially zero, so that AP core = -AP gI i,\. We can then use 
the relation 



a(f) 



Pgnd(t") ~ fgridC*" ') 
M core Af 



(B.8) 



to determine the acceleration in this time step. Finally we take 
the effects of the global acceleration of our frame into account 
in a second step, by adding 



p'(x',y',z',t) =p(x,y,z,t), 
P'(x',y',z',t) = P(x,y,z,i), 
v'(x',y',z', t) = v(x,y,z,t) - ate. 



kin 



(x',y',z',t) = ekin(x,y,z,t) ■ 



v 2 /2 + (v z - at) 2 /2, 



g'(x',y',z',t) = g(x,y,z,t)-ae z . 



(B.6) 



-a(f)At = -Av" C0 



(B.9) 



to the hydrodynamic velocity in each zone of the grid, in 
essence performing a Galilei transformation to an instanta- 
neous inertial frame in which the neutron star is again at rest. 



From relations dB.41 — dB76i . it is easy to see that the equa- 
tion of mass conservation JB. It does not change in the acceler- 
ated frame. The momentum equation in this frame is 



p'(^ + (v'-V')v'] + V'f 



■Pg 



(dv x 



8Vy 



-e r 



Note that in contrast to Eq. jB.2l > there is now an additional 
term on the right hand side, which affects the momentum com- 
ponents perpendicular to the direction of acceleration. Thus for 
instance the x-component of the time derivative of the velocity 
is 



dv' x 
~dt 



1 dx' 



^ , dv'x ^ , dv 'x 



1 dP' , dv x ,„ 

-^ + ^~ at aT (B - 7) 



where the additional (rightmost) term is negligible compared 
to v' z (dv' x /dz r ), as long as \at\ « |v z |. 



Appendix C: Explosion energy 

The explosion energy of neutrino-driven supernovae consists of 
two major contributions. The first is the recombination energy 
of the matter in the gain layer at the onset of the explosion. This 
matter consists of free nucleons and alpha particles at the time 
the explosion starts. Almost all of this mass (except for some 
fraction in the downfiows, which is accreted onto the neutron 
star) ends up in a dense shell behind the expanding shock. As 
the shock propagates outward, the temperature in this expand- 
ing shell decreases and the matter recombines to a-particles 
and later to nuclei. 

Figure IC~T1 displays the available recombination energy of 
the matter in the gain layer at the time of the explosion, 



? gam 



(^exp) 



-f 



erec(r, /exp) dV. 



(CI) 



gain layer 
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Fig. C.l. Available recombination energy, JSfJr, as a function 
of the mass in the gain layer, AM gam , at the time of explosion 
for the models of Tables IA.lHA.5l The slope of this approxi- 
mately linear relation corresponds to about 5 MeV per baryon 
(dotted line). 



Fig. C.2. Explosion energy after the recombination of the 
ejecta, E exp {t rec ), as a function of the available recombination 
energy in the gain layer at the onset of the explosion, E^(t exp ), 
for the models of Tables IA.lHA.5l For low explosion energies 
the two quantities agree well. 



Here e KC (r, t) denotes the density of recombination energy 
available when matter consists of nucleons, a-particles and 
some mass fraction of heavy nuclei, 

e rec (r, f) = B h «™ x (r, f) - (B a n a (r, t) + B h n h (r, f)) , (C.2) 

with Z?h and B a being the binding energies of a representative 
heavy nucleus (Zh, Ah = M, + Zh) from the iron group (as as- 
sumed in our equation of state) and of a-particles, respectively, 
and n™ ax = min(« p ot /Zh, n^/Nh ) and «h are the maximum and 
current number densities, respectively, of this heavy nucleus 
when «p 0t and n™ are the total (bound+free) number densities 
of protons and neutrons. 

Figure IC.ll shows that for all models of Tables IA.lHA.5l 
£rec in (?ex P ) ~ ^ ain (f exp ) x 5 MeV, when JVf in is the total number 
of baryons in the gain layer. This means that due to the partial 
assembly of free n and p in a-particles at the time of explo- 
sion, about 5 MeV (instead of > 8 MeV) remain available for 
being released by recombination during the subsequent expan- 
sion and cooling. 

This recombination is essentially complete when the shock 
has reached a radius of 3000 km (recombination to a-particles 
happens even much earlier). We denote this time by f rec . 
Fi gure IC2l demonstrates that the explosion energy at time f rec , 
£exp(frecX roughly equals the available recombination energy, 
E^. (iexp), at the onset of the explosion. This means that neu- 
trino heating essentially has the effect of lifting the total en- 
ergy of mass elements in the gain layer close to zero (i.e., 
ftot = eidn + %t + e grav ~ 0) and thus makes this matter un- 
bound and enables its expansion in the gravitational potential 
of the forming neutron star. The excess energy of this matter 
at time t rec , i.e. E exp (t lec ), is provided by the recombination of 
nucleons to a-particles and finally to iron-group nuclei. Only 
in case of higher explosion energies, E exp (t KC ) is clearly larger 
than ESrfep) (Fig. IC.2t . In this case neutrino heating in the 
gain layer is stronger and the heating time scale of the matter 
there shorter than the expansion time scale when the shock be- 



gins to accelerate outwards. Therefore neutrinos are able to de- 
posit "excess energy" in the ejecta before this matter has moved 
out of the region of strong heating. 

The second contribution to the explosion energy comes 
from the neutrino-driven baryonic wind which sets in after the 
surroundings of the nascent neutron star have been cleaned 
from the initially heated gas. Indeed this wind is an important 
energy source at "late" times. To demonstrate this, we com- 
pare the time derivative of the explosion energy, dE &xp /dt, with 
the wind power, L w i„d, and the net energy loss/gain rate L s h ck 
at the shock (Fig. IC.3i . The curve for dE exp /dt in Fig. IC.3I is 
calculated as the numerical derivative of the energy integral 

£ex P (0= f e tot (r,f)dV, (C.3) 
Jv+ 

where the integration is performed over the volume V + , in 
which the total energy e tot (r, f) is positive (see also Eq. IA.5I . 
For t > f rec this volume fills the region between an inner bound- 
ary at r ~ 200 km and the shock (except for some parts of the 
accretion downflows, where e tot may still be negative). 

The explosion energy is subject to changes by !PdV-work 
performed at, and by energy fluxes through the boundaries of 
V + , in particular by the wind, whose power is given by the sur- 
face integral 

^wind — Q) (etot + free + 9) max(v r , 0) dS. (C.4) 

Jr=200km 

This expression takes into account the total energy (e tot = pe tot 
with e tot defined by Eq. IA.6i of the wind material streaming 
through the inner boundary radius into V + , the energy that will 
be set free at larger radii by recombination (Eq. IC.U . as well 
as the work performed by pressure forces. Here we have ne- 
glected effects due to downflows by omitting contributions to 
the surface integral from zones with negative radial velocity. 

The change of the explosion energy due to energy flow 
through the outer boundary of V + (i.e., the shock), is given by 
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Fig. C.3. Evolution of the time derivative of the explosion en- 
ergy (dEexp/df, thick solid) for Model B 18. Also shown are the 
wind power at a radius of 200km (L w i nc j, dotted), the energy 
loss/gain rate at the shock by fdV work and swept-up mat- 
ter (L s hock, dashed), and the sum of the latter two quantities 
(Avind+shock, thin solid). L wind+shock agrees well with dE exp /dt 
for t > f rec (right of the vertical line). 

the net energy loss/gain rate 

^shock = (I) [(^tot + free 

+ T)v r +(e tot + e rec )fl s ]dS. (C.5) 

Jr=RJfi)+ 

The integration has to be performed over a surface located 
slightly upstream of the shock. Compared to Eq. ( IC.41 an ad- 
ditional term arises here from the motion of the shock, which 
propagates with a local velocity R s (8). 

Figure lC31 shows that these two terms explain the evolution 
of dEexp/df for t > f rec , i.e. 

dEexp/df « L w i n d + L s hock (C.6) 

holds at late times, and the thin and thick solid lines in Fig. IC.3l 
almost coincide. Note also that L w i n d s> |L s hockl- This is true 
for all our models, and therefore the increase of the explosion 
energy after about 0.5 s post bounce is (almost exclusively) as- 
sociated with the time-integrated wind power (see Fig. lC.4> . 

The relative importance of the two major constituents of the 
explosion energy that we have discussed here, i.e., the nuclear 
recombination energy of the matter in the gain layer and the in- 
tegrated power of the neutrino-driven wind, varies with the ex- 
plosion energy. In our "standard boundary contraction" models 
the fraction of the explosion energy provided by recombina- 
tion drops from about 70% for the low-energy models to about 
30% for the model with £ exp (l s) ~ 1.5 x 10 51 erg (Fig. PI. 
This fraction declines because the wind power is proportional 
to a higher power of t he luminosity (L w ind « L" with a » 3; 
iThompson et al.l200lh and although the mass in the gain layer 
at the onset of the explosion scales linearly with the boundary 
luminosity (Fig. llOt . 

For the "rapid boundary contraction" cases the wind con- 
tribution is even more important, e.g. for Model W12F-C 
£rec in ( f exp)/£exp(l s) ~ 0.2, i.e. about 80% of the explo- 
sion energy are generated by the neutrino-driven wind. For a 



Fig. C.4. Relation between the increase of the explosion en- 
ergy between t = 0.5 s and t — 1 s, AE'^ 5s > an d the integrated 
wind power during this time interval, A£fj?'? s , for the models 
of Tables |A~THA31 



1.0 




0.0 L , , , , . , , , , , , , , , 
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Fig. C.5. Ratio of the recombination contribution to the total 
explosion energy 1 s after core bounce, £rec m ( f exp)/£exp(l s), as 
a function of E exv {\ s) for the models of Tables lAJTlA. 51 For 
low-energy models the recombination contribution dominates, 
whereas for higher explosion energies the wind contribution 
becomes more important. 

fixed boundary luminosity the wind power is higher in this 
case than for the "standard boundary contraction", because 
Ewind increases with d ecreasing neutron star radius (see e.g. 
Thompso n et alJl200ll) . However, AM ga i n (f eX pX and thus also 
^fec'VexpX are similar for models with "standard" and "rapid" 
boundary contraction and the same L^. This is so because two 
effects compensate each other roughly: On the one hand the 
density at a given radius r in the gain layer is lower for a faster 
contraction (p r (r, f exp ) < p s (r, fg Xp ), where r and s denote the 
rapid and standard contraction cases, respectively), but on the 
other hand also the gain radius is smaller and thus located in a 
region of higher density, p r (R r g , f exp ) > p s (^ g , f exp )- 

Figure IC6l indicates that the explosion energy is still in- 
creasing at t — Is when we stopped most of our simulations. 
Yet, with the subsequent drop of the core luminosity (we as- 
sume a r 3/2 behaviour at t > t L , see Eq. ID.11> also the wind 
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time [s] 

Fig. C.6. Evolution of the explosion energy for Models R18- 
c and W18-C which are listed in Tabled The rotating model 
R18-C attains an explosion energy which is higher than in the 
non-rotating case W18-C due to a larger gain layer mass at the 
time of explosion. Note also that the explosion energy in both 
cases is still increasing at t — Is. 

power, which is proportional to L" (see above), must decline 
strongly. Therefore the explosion energy will grow only mod- 
erately. In case of the long-time simulation B18-K it rose from 
1.14xl0 51 ergatf = 1 s to 1.43xl0 51 erg at t = 2 s, and reached 
1 .46 x 10 51 erg by the end of the simulation at t — 3.6 s. 

Appendix D: Neutrino transport 

D.1. Transport equation 

We start from the equation of radiation transport in spherical 
symmetry 

1 ° I+ ° I+ l _Zllp = S> (D.l) 
c ot or r op 

where / = I(t,r,e,p) is the specific intensity, S = S(t,r,e,p) 
is the source function, e is the neutrino energy, p = cos and 
9 is the angle between radiation propagation and radial direc- 
tion. Solid angle integration yields the zeroth angular moment 
equation, 

13 1 d , m 1 f +1 

c ot r L or 2 J_j 

with {J, H}(t, r, e) := | J^ 1 dp p^°'^I(t, r, e,p). Integration over 
energy leads to 

-E+-—(r 2 F) = Q + -Q- (D.3) 

ot r L or 

with {E, F}(t, r) :— 4n de {J/c, H}(t, r, e) being energy den- 
sity and energy flux, respectively. The source term has been 
split in an emission rate Q + and an absorption rate Q = K. d cE, 
which is proportional to the energy density. The flux factor is 
defined as the ratio of flux to energy density, 

f(r,t) := F(r,t) / cE(r,t). (DA) 



In neutrino transport simulat ions solv i ng the full 
Boltzmann equation (see e.g. iBuras et afll2003t l2006albl) this 
quantity shows only little short-time variability during most 
phases of the supernova evolution. Therefore df/dt = is an 
acceptably good approximation. With L :- 4nr 2 F = 4nr 2 fcE 
one can now rewrite Eq. iD.3l as 

-L + c eff -L = 4nr 2 c eS {Q + -Q-}, (D.5) 
ot or 

where an effective speed of neutrino propagation has been in- 
troduced as c e ff := cf. Provided c e ff were known, the solu- 
tion of Eq. iD.5\ requires considerably less effort than the nu- 
merical integration of Eq. JD. li . For vanishing source terms 
Q + and Q~ the neutrino energy or number density is just ad- 
vected along characteristics r(t) — ro + c^t. Although c e ff 
depends through f(r, t) on the solution of the transport prob- 
lem (Eq. ID.4> . neutrino transport calculations in the neutrino- 
decoupling layer of forming neutron stars reveal that it can 
be well fitted by a r-dependent functio n which depends o n 
the steepness of the density profile (see Ijankalll990t Il991bl) . 
Assuming further that the (medium-dependent) coefficients Q + 
and k = K.Jf - 4nr 2 Q^/L are constant between two points 
(r, t) and (r*, £*), which are connected by a characteristic line, 
i.e., 

r* =r-c eff (t-f*), (D.6) 
Eq. (ID.5> can be integrated analytically to yield 

L(r,t) = L(r*,t*)e-~ KC ^ if -'* ) 

+^®L ([i _ e-fc-rCW*)] U + ( ~ Kr * _ 1} 2] 

+Kc eS (t - t*)[2kr* + A-c eff (f - f*) - 2]}, (D.7) 

where L(r, t) and L{r*,t*) are the luminosity values at both 
ends of the characteristic line. 

We use Eq. iD.7\ to construct a numerical scheme to solve 
Eq. iD.5\ in the general case: We assume that the luminosity 
is known at the cell interfaces of a one-dimensional radial grid 
for a time f"~', and that the cell-averaged values of the quan- 
tities needed to compute the emission rate Q + and absorption 
coefficient k are also known for that time. As a further simplifi- 
cation we do not allow neutrinos to propagate in negative radial 
direction (actually this is granted by defining a non-negative 
function for the flux factor, see Sect. ID. 3\ . Then the luminosi- 
ties at t" = t"~ l + At for each zone interface (starting at the 
innermost zone) can be computed using Eq. (ID.7> . In doing 
so we have to distinguish between two cases (see Fig. lD.U : If 
c e tf At > Ar, we can use point A as the starting point of the 
integration, (r*,t*) = (r,_i,fA). The luminosity at this point 
is derived from a linear interpolation between L(r,-_i, f"~') and 
L(r,_i, f") (which is already known, as we are integrating out- 
wards). If c e ffAf < Ar, we use point B, the luminosity at this 
point being given by a linear interpolation between L(r,_i , f"~') 
and L(r,-, ?"-'). 

For time integration we use a predictor-corrector method: 
The transport routine is called two times. In the first (predic- 
tor) step the luminosities, emission rates and absorption co- 
efficients of the last time step [L" -1 , Q n_1 , k"' 1 ] are used to 
compute preliminary values {Q",k") for the neutrino-medium 
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Fig. D.l. The solution at (r',t") is computed from the data 
at a point (r*,t*) located on the same characteristic line. 
Depending on the grid spacing, Ar, the time step, At, and the 
effective speed of neutrino propagation, c e g, either point A or 
point B must be used. The solution there can be obtained by 
interpolation in time or space, respectively. 

coupling at the next time level. In the second (corrector) step 
the final values [L n , Q n , k"] are calculated using [L n ~ l , \(Q"~ l + 
Q n ), + £")] as input. 

Equation ( ID.5> is solved not only for the energy luminosity 
L — L e , but also for the number luminosity L„ = Anr 2 F n = 
4nr 2 fcn (n is the particle density and / is assumed to be the 
same flux factor as for the energy transport). Furthermore the 
equation has to be integrated for three neutrino types, v e , v e , 
and v x (the latter denoting v^, v^, v T , and v T , which are treated 
identically). In the following we will use indices v e {v e , v e , v*} 
and a e {e, n) to distinguish between these different cases. 

In the 2D case the neutrino transport is assumed to proceed 
only radially, i.e. lateral components of the neutrino flux are 
ignored and Eq. flD.5l > is integrated independently on different 
radial "rays", i.e. in radial direction for every lateral zone of the 
polar coordinate grid. Total luminosities of the star are obtained 
by summing up the flux densities L/Anr 2 for all angular cells 
(at a given radius r), appropriately weighting them with the 
corresponding surface elements. 

D.2. Boundary conditions 

To integrate Eq. JD.7> outwards, time-dependent boundary con- 
ditions are required for the luminosities L e v . and L„ v ., where 
v; = v e , v e , v x . We assume L e v . to be constant for a time interval 
ti (typically 1 s), and to decay subsequently with a power-law 
dependence in time: 



where 



hit) 




L^°K Vc h(t), 
L^°K % h(t), 
Lf'° K v h{t), 



if t < t L , 
if t > t L . 



(D.8) 
(D.9) 
(D.10) 

(D.ll) 



The constants K v . denote the fractional contributions of the in- 
dividual luminosities to the total neutrino luminosity. They ful- 
fill the requirement 
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Fig. D.2. Evolution after core bounce of the neutrino luminosi- 
ties at a Lagrangian mass shell of 1 . 1 M fro m a supernova sim - 
ulation with Boltzmann neutrino transport jBuras et"al1 Eb03). 
After an initial phase of 50 ms duration, the sum of the v e 
and v e luminosities as well as the v^jvj luminosities vary only 
slowly. 



The functional form used in Eq. ( ID. 1 11 can be motivated 
by the Boltzmann transport calculations of Bura s et all J2003). 
These show that after a transient phase of ~ 50 ms, which is 
short compared to the explosion time scales of our simulations, 
the sum of all luminosities is almost constant or varies only 
very weakly, at least over the next ~ 250 ms, for which data 
from the Boltzmann transport simulations are available (see 
Fig.lD~2l. 

According to Eqs. ( ID.8HD.rTl we need to prescribe the 
time scale ti and the total initial luminosity l}° x,a . However, 
instead of choosing these two quantities as basic parameters 
of our models, we prefer to prescribe ti and the gravitational 
binding energy AE™ col . e that is released by the neutron star 
core asymptotically (i.e. for t — * oo) via neutrino emission. 
Introducing the energy that the core looses up to time t 



A^y.coreW 



Jo 



(D.13) 



the following relations hold for the asymptotic energy loss 



ae: 



Jo 



4°'' /!(f)df : 



3A2?" (t L ) 



Ky+K, +4Ky=L 



(D.12) 



(D.14) 

i.e. our ansatz of Eq. ( ID. 1 1> implies that 1/3 of AE£° core is ra- 
diated away within the chosen time interval ti in neutrinos and 
antineutrinos of all flavours. 

We also prescribe the mean energies of neutrinos entering 
the computational grid at the inner boundary. The correspond- 
ing values are chosen to be <e Vc ) ib = 12MeV, <ev c ) ib = 16MeV, 
and (e Vi ) lb = 20MeV, and kept constant during our simu- 
lations. Thereby also the number fluxes L„ v . = L e>Vi /(e Vi ) at 
r = are defined. 

The total lepton number lost by the neutron star core un- 
til time f, normalised to the total baryon number A%, cole of the 
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core, is given by 



with the Fermi integrals defined by 



AFe,core(0 = N b \ ole f (L n ^(R ih , ?') - L„^(R ih , t')) At'. 
Jo 



(D.15) 



T tot,0 . / ^ 



(D.16) 



For t — ?l this yields 

AF e , core fe)= Nbcoie \ {€v ^ {e ^ j 

We assume that the lepton number loss during time interval ?l 
is proportional to the energy loss during this time. Therefore we 
choose K Y . = const, because (e Vl ) lb = const, and set K Vc = 0.2, 
K Ve — 0.215 for the calculations in this paper. K Vx follows from 
Eq. flXHt , 

D.3. Neutrino distribution function 

To calculate the source terms in Eq. (ID.5> or JD.7> we have 
to make an assumption about the neutrino energy spectrum, i.e. 
about the energy dependency of the specific intensity, which for 
particle energy and particle number is linked with the particle 
distribution function v in the following way: 

J2,3) ' 



Iy,[„, e }(t,r,e,fi) 



(he) 3 



c fr, v (t,r,e,fr), 



(DAI) 



where the exponent of 2 applies for number transport and the 
exponent of 3 for energy transport, corresponding to the indices 
n and e, respectively, of I v . We assume that /b,v can be written 
as product of a Fermi-Dirac distribution function, 

^ )= i + «U-0 - (D - 18) 

and an angle-dependent function g v , 



/ Av (r, t, e,jX) = g v (r, t,fr) f FD 



ij Y (r,t)\, (D.19) 



\k B T v (r,t) 

where in general the spectral temperature and degeneracy pa- 
rameter, T v and J] v , are different from the matter temperature 
and equilibrium degeneracy parameter. 

Furthermore we assume that r\ v is just a function of the op- 
tical depth r„: 

riv(T v ) = 77eq, v (1 - e _Ty ) + rjo.y e~ T ", (D.20) 

where rf e ^ v is the equilibrium degeneracy parameter and r]Q y y is 
a chosen value as typically found in detailed transport calcula- 
tions for Ty — > 0. The values for the different neutrino types are 
(cf.lJankall991bHjanka & Hillebrandll98llMvra & Burrows! 
1990; lKeilet all2003l) : 

^eq ; v t = (y"e + j"p ~ f-l n )/k B T, TJq^ = 3, 

^eq.v. = -?7eq,v c , ?70,v c = 2 , (D.21) 

?7eq ; v, = 0, 770,y t = 0. 

With j] y defined, T v can now be calculated from the local 
average neutrino energy, which is computed from L e v and L„ jV 
as 

(e„) = L ejV /L„ jV = Ey(r,t)/n v (r,t) 
r+l 



poo p+ 1 

J de J j d/dl„ t y(r,t, e >y u) 
= k B T v r3(i7y)/r 2 (r]y) 



(D.22) 



T n (rf)= I dx^/ FD (x,/7) 
Jo 



(D.23) 



Thus the energy-dependent part of f D v is fully defined. The 
angle-dependent part is related to the flux factor by 



fy = Fy/cEy = 



f-1 d ^^ X^'dyUgy^) 



= ( f ly). (D.24) 



To solve Eq. iD.5i only f v is needed, not the angle-dependent 
function g v (fr). Far outside of the neutrinosphere, f v is ap- 
proaching the vacuum solution. The latter can be derived under 
the assumption that neutrinos are emitted isotropically from the 
sharp surface of a sphere with radius R v , which is located at a 
distance r from the observer. In this case the flux factor is 



fv,\ 



[l+ Vl-CRv/r) 2 ]. 
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/v,vac approaches 1 for r — > oo (free streaming limit) and 
/v.vac = 1 /2 at the neutrinosphere. In a more realistic situation 
the neutrinosphere is not a sharp surface but a layer with finite 
thickness in which neutrinos gradually decouple from the stel- 
lar medium. In detailed transport calculat ions f v (R v ) is there- 
fore found t o be about 1/4. (see e.g. lJanka & Hillebrandll98^: 
Janka Il990h . How fast f v approaches / v , va c (with declining op- 
tical depth) depend s on the steep ness of the density gradient at 
the neutrinosphere ll.Tankall990h . Inside the neutrinosphere de- 
tailed transport calculations show that the flux factor behaves 
roughly like f y oc t 7 " with m < 0. 

Taking all this into account, the following function consti- 
tutes a good approximation for the flux f actors from detailed 
transport calculations (Janka 1991b, 1990): 



fy(Ty) = 



ki+D] 



1 +(1 +D)(1 

i/4(T V /T„ir, 



d y' 



if 
if 



Ty < r Vj i, 



> T v ,i, 

' (D.26) 

Here D = yj 1 — (R v / r) 2 , the neutrinosphere radius is defined 
by Ty(Ry) = r Vj i and we adopt t Vj i = 1.1. The power-law index 
m is chosen such that / v (10) = 1/25, and n is defined by a local 
power-law fit of the density profile around the neutrinosphere, 
p(r) oc r~ n . A higher value of n therefore means a steeper den- 
sity gradient. 

DA. Neutrino reactions 

For calculating the neutrino-matter interaction rates the follow- 
ing reactions are taken into account: charged-current processes 
with neutrons (n) and protons (p), 



v e + n =i p + e 



(D.27) 
(D.28) 



thermal electron-positron (e*) pair creation and annihilation, 
e + + e~ =i v, + v ; (i = e,//, t), (D.29) 
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and neutrino scattering off nuclei (A), nucleons, and electrons 
and positrons, 



v,- + A ^ Vj + A, 





In] 






H 










IPJ 







v; + e ^ v,' + e 



(D.30) 
(D.31) 
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D.5. Optical depth 

Knowledge of the optical depth is necessary to evaluate 
Eqs. (ID.20> and (ID.26> . For this purpose it is sufficient to com- 
pute t y approximately by considering only the most relevant 
neutrino processes and assuming that the neutrino spectrum is 
given by the spectrum for local thermodynamic equilibrium. 
This means that instead of Eq. jD.19t we use 



■O e v, r) = /fd 



k B T(r) 



?7eq,v0) 



(D.33) 



with 77eq. v and T instead of tj v and T v . 

The "transport optical depth" is defined as the integral 



TtyO) 



-I dr ' 



<*tvXO 



(D.34) 



of the energy-averaged "transport opacity" (i.e. the opacity 
which is relevan t for momentum transfer), (K t . v )(r) (see, e.g., 
IStraumamill 19891 burrows & Thompsonll20 04). In the follow- 
ing, all neutrino interactions included in evaluating the opacity 
are calculated without final-state lepton blocking, unless other- 
wise stated. 

The most important opacity-producing reactions are scat- 
tering off nucleons (n,p) and nuclei (Zj,Aj), where j — 1,2,... 
denotes the considered nuclear species, and absorption by neu- 
trons and protons in case of v e and v e , respectively. Thus one 
can write 

<*,,v> = <#£>+ J] <<•;>■ (D.35) 

ie{n,p3j) 

Here the (neutrino-flavour independent) scattering opacities are 
to lowest order in neutrino energy over nucleon rest mass (i.e., 
without effects of nucleon recoil, thermal motions, and weak- 
magnetism corrections): 



«; p > = 



-or 2 + (Cy - l) 2 



cro 



(m e c 2 ) 2 



<e v > « p 



5a 2 + 1 



0"o 



24 (m e c 2 ) 2 
O = \ A )\- C a -1 + t-P-Ca- C v )f 

X (^¥^ > ^' 



(e v ) n n , 



(D.36) 
(D.37) 

(D.38) 



for scattering off protons, neutrons, and nucl ei with number 
densi t ies n p , n„, and n A , ^respectively (see, e.g..|Freedman et all 
ll977HStraumannlll989tlBurrows & Thompsonl2004 . The ab- 

sorption opacities for v e and v e are 



<4> 



(<e 2 >+2A<e Vc > + A 2 ) 0«e Ve », 



i 2 °~0 

4 i3a + 1} (^)2 "p 



(<6?>*+2A<ev c >* + A 2 <e°>*). 



(D.39) 



(D.40) 



Tsee lTubbs & Schrammll9"75llBruennl 19851) . 

Here <r = 4G\m 2 ri 2 Inc 2 = 1.76 x l(T 44 cm 2 (with the 
Fermi coupling constant Gf), A = 1 .2935 MeV is the rest mass 
difference of neutrons and protons, a = 1.254, Ca - \,Cy - 
i + 2 sin 2 6 W , and sin 2 % = 0.23. 



In deriving Eqs. JD.36t - JD.4CH (as well as for all rates 
and source terms given below) the electron and positron rest 
masses are ignored (m e c 2 <k e v ) and nucleons and nuclei are 
assumed to have infinite rest masses (m n ^.c 2 » e v ) and to be 
nondegenerate. For electrons, phase space blocking is included 
in Eq. JD.40> by the factor 



/ (e v ) + A 
knT 



(D.41) 



which accounts for the fact that a significant fraction of the pos- 
sible final electron states may be occupied. Phase space block- 
ing can be neglected in k? (Eq. ID.40> . because the positrons 
are non-degenerate. 

The neutrino energy moments are (generalising Eq. ID.22> 
given by 



«> = (k B Ty) 



«>* = (k B T y ) 



r 2 (r] v ) 

„ T 2+ „ (ilv-A/k B T v ) 

r 2 (T]y) 



(D.42) 
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and for evaluating Eqs. ( ID.36> - db.4m to compute r tv 
(Ea. ID.34l for use in Eq. (ID.26> . we take tj v = j] &q v and T v = T. 



In contrast, Eq. JD.20> is evaluated with the "effective opti- 
cal depth for equilibration", 



dr'<K effjV )(r'), 
where the effective opacity is defined as 

<Keff,y) 



+ <y>- 



(D.44) 



(D.45) 



Here the spectrally averaged absorption opacity, (rf), is taken 
to include neutrino-pair annihilation to e ± -pairs (Eq. ID.29> . 
which is assumed to be the most important reaction for pro- 
ducing v x v x pairs. Both (^) and (rf + k^ v ) are evaluated for the 
"true" (not the local equilibrium) neutrino spectrum (i.e. for the 
spectral temperature T v and the spectral degeneracy j] v instead 
of T and ?7eq,v) by employing the source terms from the neutrino 
transport solution of the last time step. 
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D.6. Source terms 
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These rates hold for neutrinos v or antineutrinos v of all 



Solving Eq. iD.H requires the knowledge of the emission rates, 
Qy., and absorption coefficients, k v = i$lf v , which appear in 
this equation. Since Eq. (ID.7> is used to determine the number 
fluxes, L„ zV , and luminosities, L e >v , of all neutrinos and antineu- 
trinos v € {ve,v e ,vj, the source terms need to be calculated 
for the neutrino number as well as energy. In the following, all 
these neutrino source terms are derived without taking into ac- 
count final-state lepton blocking, unless otherwise stated. As in 
Eq. JD.45I I. k v is defined to include the contributions from the 
/^-processes, Eqs. JD.27> and JD.28> for v e and v e as well as 
those of e + e~ pair annihilation (Eg. ID. 291 . The absorption co- 
efficient can be computed from the corresponding neutrino 
absorption rate by 

< = Q-A7rr 2 f v IL v = (Ql + ) 4nr 2 f v /L v . (D.46) 

For the number transport the neutrino absorption and emis- 
sion rates (in units of number per cm 3 per second) by charged- 
current /^-reactions between leptons and nucleons can be writ- 
ten with our approximations for the neutrino distribution func- 
tion and the appropriate statistical weights for the leptons as 
follows: 



V e 
ft? 

K 



L ev n R <e2 ) + 2A(e Vc ) + A 2 



= ac ^K 0(<O) ' 

L e v c n p {el)* +2A(e v X + A 2 (e°>* 
= <TC : 

4nr z cf Vc 



= -crcn p n e - + 2A<e e ->* + A 2 <4>*], 

= icrcn n n e+ [<4) + 2A(e e+ ) + A 2 ], 
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where cr = j(3a 2 + l)o"o / (m e c 2 ) 2 and the electron (positron) 
number density is 

8tt 



(he)- 



The electron and positron energy moments are given by 

,„ T2+«(?7e) 



<<£) =(k B T)' 
<e e ">* =(k B T)« 



r 2+n (rj e -A/k B T) 



(D.51) 

(D.52) 
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The annihilation and production rates of neutrino num ber in 
e + e~ pair react ions are given by (adapted fromlSchinder et all 
Il987l see also |jankal[l991al and Ijankall 1 99 lbl and references 
therein): 



o-qc 



prod 



(4;rr 2 C ) 2 <e v ><e„) 

1 1 
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6 

1 



2<t>(f V , Xv )C 2 A y+Cly 



2\2 



/v/> 
cr c 



fvfv (is-.! 



18 (m e c 2 ) 2 



n e -«e 



+ -(m e c 2 ) 2 (2C 2 v - 



9 fyf, (m e c 2 ) 



ci v ) 



(Cly + C 2 )(e e -)(e e+ ) 



C Av ) 



flavours. In Eq. (ID.46> . Hy and 7^™ have to be used instead of 
Q v and Q y nn for computing the absorption coefficient for the 
number transport. In Eq. jD.54l i. €>(/ v ,^ v ) is a geometrical fac- 
tor, 



1 - 2/ v /v +XvXv + ^(1 -Xv)( l ~Xv) 



(D.56) 

where we express the variable Eddington factor Xv in terms of 
the flux fa ctor f v (E q. ID.4> using a statistical form, which was 
derived by MinerboJ (119781) on grounds of maximum entropy 
considerations (for photons or nondegenerate neutrinos, as as- 
sumed here): 

0.01932/;, + 0.2694 f 2 



Xv = <Mv) 



1 

3 + l 



■ 0.5953 f v + 0.02625 f 2 ' 



(D.57) 



The weak coupling constants in Eqs. jD.54l and ( ID.55I > are 

given by 



for v e 
for v e 



Cv 



+5 + 2sin 2 6>v 
-i + 2 sin 2 # 



'w 



for v e 
for v E 



IVe, V e 
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The source term which describes the rate of change per unit 
of volume in the evolution equation of the electron lepton num- 
ber of the stellar medium is 

N = y e "b = (K - K) - (K - *?.)• ( D - 6 °) 

The source terms which account for the absorption and 
emission of energy through v e and v e are computed in analogy 



(D.50) to Eqs. dD~?7l - dD30l as 



Ql 



- crc 



L e , Vc n n <e 3 ) + 2A<6 2 ) + A 2 (6 Vc ) 
Anr 2 cf Vc (e Vt ) 



■0«e Vc », (D.61) 



= crc 



Anr 2 cf Vc 

<e 3 )* + 3A(e 2 )* + 3A 2 (e v X + A 3 (e° )* 



Ql = y «p» e - [<e e 3 ->* + 2A<e 2 _)* + A 2 <e e ->*] , 

Ql = ^»n»e + [(4> + 3A(6e 2 t )+3A 2 (6e + )+A 3 ], 
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The annihilation or production of energy in neutrinos (v) by 



e + e pair reactions is given as ( Janka 1991a) 



Ql 



cr () c 



2^(fy, X y)C 2 A y+Cly^ : 
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■cL) 
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Table D.l. Weak coupling constants for v and v scattering off 
e + or e~ (cf. Ea. lrT67l . C* stands for C* = (C A - 1 ) 2 - (C v - 1 ) 2 , 
Ca = j, Cy = 5 + 2 sin 2 6w, and can be or v T . 
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c 3 
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(C v - C A ) 2 


c 2 -c 2 
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(C v - C A ) 2 


(C v + C A - 2) 2 




Kv e 


(C v - C A ) 2 


(C v + C A - 2) 2 




y v e + 


(C v + C A - 2) 2 


(C v - C A ) 2 


c* 



For annihilation of antineutrino (v) energy, (e 2 ) has to be 
replaced by (e?) and (e v ) has to be exchanged with (e v ) in 
Eq. (ID.651 . while the production of v and v was assumed to 
be symmetric and both rates are given by Eq. JD.66> . 

Also in scattering processes energy can be exchanged be- 
tween neutrinos and the stellar medium. F or scat tering off e~ or 
e + , using the rates of Tubbs & Schramm ( 1975), and ignoring 
electron phase space blocking in the final reaction channels, the 
following spectrally averaged expression for t he energy tr ans- 
fer rate per unit of volume can be derived (see Janka 1991b): 

n - —(c j. ir Lg - V 
Uvt ~ l2 (Ll 6 2 VeC 2 ) 2 " e 47n-2 C /v<fv> 



[<e 2 )«6e> 



+ -^m e c ) 



3 

H — 



(m e c ) [(£ v ) - — 



c, + he 
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where e can be e + or e~ and v stands for neutrinos or antineutri- 
nos of all flavours and the constants C\ , C2, C3 for the different 
combinations are listed in Table ID. ll The term 3ra e c 2 /4 in the 
bracket results from a merge of the rate expressions for the lim- 
its of relativistic and non-relativistic electrons. In the latter case 
the neutrino-electron scatteri ng cross sect ion is proportional to 
e v /(m e c 2 ) for e v » m e c 2 (cf. Se hgalll974l) . 

Every transfer by neutrino-nucleon scattering, which is 
only "nearly conservative", is taken into account following 
lTubbsldl979t) . The corresponding rate is (see Janka 1991b): 



QvN =T 



1 o-qc 



4 (m e c 2 ) 2 



4nr 2 cf v {e v ) 
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with 



n P _ j |[(C V - l) 2 + § a 2 ] forN = p, 
U(l+5ar) forN = n. 

The symbol v stands again for neutrinos and antineutrinos of 
all flavours. Also scattering contributions are included in the 
energy generation rate Q + and absorption coefficient k used 
in Eq. jD.7> . Considering scattering as an absorption process 



followed immediately by an emission process, we add the net 
energy exchange rates Q ve - , Q ve * , Qvp an d Qva to Q~ (used for 
computing k v in Eq.[^2} when the rates are positive (i.e. in case 
of energy transfer from neutrinos to the stellar gas), and the ab- 
solute values of Q ve -, Q ve +, Q vp and Q yn to otherwise. The 
total neutrino energy source term to be used in the gas energy 
equation including the contributions from v e and v e absorption 
and emission, vv pair production, and all scattering reactions is 

Qe = X (qi-qi)+ X (er-C d ) 

ve (v e ,v e ) ve {v e ,v^,v T } 

+ £ {Ql e+ + Ql e - + Ql P + Q s yn ), (D.69) 

where = + and Q v f = (%* + Qf d . 

In practise, however, the lepton number source term 2n as 
well as the energy source term for the hydrodynamics part of 
the code is not computed from Eq. jD.60t and jD.69> . respec- 
tively, but from the luminosity change between points (r',f) 
and (r*,f*) (cf. Fig. lD.H . The source terms Q' N and Q' E for a 
grid cell i at time level f are then given by 



L™(r\f)-L^(r*,t*) 
AV~ 



L l °\r\t n )-Lf{r*,t*) 
AV> ' 
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where AV, = y (rr - r ) is the part of the cell volume crossed 
by the characteristic line between (/, t") and (r*, f*), L}° 1 is the 
sum of the luminosities of neutrinos and antineutrinos of all 
flavours, and L dlff is the difference between the v e and v e num- 
ber fluxes, L v .,„ ■ 



■ L Ve7 „. Equations iD.lOh and ( ID.71> work well 
as a description of the neutrino sources in the gas equations 
only, if the neutrino fluxes do not exhibit a large degree of vari- 
ability on the radial and temporal scales of the r-t cells. This, 
however, is reasonably well fulfilled in the context considered 
in this paper. 

Finally, the outgoing neutrino fluxes transfer also momen- 
tum to the stellar fluid. To account for this, we include a mo- 
mentum source term Qm which enters the Euler equation of 
the hydrodynamics solver. It is sufficient to include only the 
most important reactions, by which neutrinos transfer momen- 
tum, i.e. v e and v e absorption on n and p, respectively, and 
the scattering processes of v and v of all flavours off nucleons 
and nuclei (pair processes and electron/positron scattering can 
be safely ignored). For a neutrino or antineutrino v, the corre- 
sponding rate (in units of erg/cm 4 ) is 



Ql 



4nr 2 c 



+ z 

ie{p,n,A,} 



<fv> 
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where the first term in the sum is relevant only for v e and v e . 
The energy averages of the scattering transport opacities, 
and of the absorption opacities, ^, all weighted by the neutrino 
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energy, are given by 
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x (<e 3 >* + 3A(e 2 >* + 3A 2 (e, c >* + A 3 (e° >*) . 

The energy moments (e") and (e")* are given in Eqs. JD.42t 
and JD.43I I. They are calculated using the nonequilibrium neu- 
trino spectral parameters T v and rj v . The momentum source 
term in the equation of gas motion then reads 



V€{v e ,V^,V T ,V e ,V^,V T ) 



(D.78) 



It was not included in the simulations presented in this paper, 
but will be taken into account in future calculations. 

The implementation of the source terms Q^, Qe, and Qm 
into the framework of our PPM hydrodynamics c ode was dis- 
cussed in detail by iRampp & .Tankal J2002I) and iBuras et a"fl 
J2006al). 

We finish by pointing out that the approximative neutrino 
transport scheme developed here employs two basic assump- 
tions, which are radical simplifications of the true situation: 

1. In deriving Eq. iD.ll from the transport equation we as- 
sumed that the flux factor f(r, t) is a known function, al- 
though it is actually dependent on the solution of the trans- 
port problem (see Ea. ID.4t . Equation iD.l\ certainly has the 
advantage of analytic simplicity, but also has a severe dis- 
advantage: The source terms can be very large and the nu- 
merical use requires a very fine grid zoning at high optical 
depths. The cell size should fulfill the constraint that the op- 
tical depth of the cell stays around unity or less. Moreover, 
the implementation of the source terms in iD.7\ and the 
medium sources (Eqs. ID. 701 ID.7U is not symmetric and 
the numerical scheme does not strictly conserve the total 
lepton number and total energy of neutrinos plus gas. 

2. For treating the spectral dependence, we made the assump- 
tion that the neutrino phase space distribution function can 
be factorised into a product of an angle-dependent function 
g v and an energy-dependent term, which we assume to be 
of Fermi-Dirac shape. This certainly constrains the spectral 
shape, but the factorisation also implies that the flux-factor 
is assumed not to be an energy-dependent quantity. This 
in turn means that the mean energy of the neutrinos flux, 
(fy)flux = L etV (r,t)/L„ tV (r,t) is identical with the mean en- 
ergy of the local neutrino density, (e v )i cai = E y (r, t)/n v (r, t). 



This is certainly a problematic simplification in view of the 
fact that the neutrino interactions with the stellar medium 
are strongly energy-dependent. 

Nevertheless, the described neutrino transport treatment 
represents a practical approximation which is able to reproduce 
basic features of more detailed transport solutions and yields 
agreement with those even beyond the purely qualitative level. 



